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ABSTRACT 


Classically,  the  T-matrix  method  is  a  procedure  to  exactly  compute  the  multiple  scattering  of  an 
incident  wave  from  a  “cloud”  of  objects,  given  knowledge  of  the  free-held  scattering  properties 
of  a  single  object  for  an  arbitrary  incident  wave.  For  acoustic  waves,  Profs.  Baker  and  Scan- 
drett  have  extended  the  T-matrix  method  to  the  case  in  which  the  radiation  sources  are  also  the 
scatterers,  that  is,  to  the  case  of  an  array  of  active  transducers.  This  thesis  is  the  first  successful 
practical  demonstration  of  the  T-matrix  method  applied  to  an  active  sonar  array  for  which  a 
finite-element  model  was  employed  to  compute  the  scattering  properties  of  a  single  transducer. 
For  validation,  a  T-matrix  model  of  a  linear  array  of  piezoelectric  spherical  thin-shell  transduc¬ 
ers  was  modeled,  for  which  analytical  approximate  values  of  the  T-matrix  element  values  are 
known.  Subsequently,  a  T-matrix  model  of  a  linear  array  of  piezoelectric  class  V  flextensional 
“ring-shell”  transducers  was  modeled.  Beam  patterns  of  the  linear  array  models  computed  with 
the  T-matrix  method  are  compared  with  those  of  an  array  of  point  sources,  demonstrating  that 
the  T-matrix  method  produces  more  realistic  beam  patterns,  especially  for  end  fire  arrays. 
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CHAPTER  1: 

Introduction 


1.1  Motivation 

In  many  sonar  applications,  it  is  desirable  to  use  a  high  power  and  directional  sonar  transducer. 
Their  use  reduces  reverberation  and  can  propagate  further  distances.  Also,  they  provide  more 
accurate  estimate  of  bearing.  These  transducers  are  used  in  arrays  in  order  to  maximize  the 
directionality  and  steer  optimum  acoustic  beams.  Figure  1.1  shows  two  examples  of  volumetric 
sonar  array  applications. 


(a)  (b) 

Figure  1.1:  (a)  A  bi-planar  array  containing  of  24  Sparton  of  Canada  600  Hz  class  V  flex- 
tensional  transducers  [1]  (b)  Bicycle  Array  composed  two  of  the  610  Hz  class  V  flextensional 
transducers  [2] 

The  interaction  between  array  elements  can  affect  array  performance.  Especially,  this  effect 
becomes  severe  for  high  power  and  directional  array  elements  because  they  are  commonly 
operated  at  a  resonance  frequency  and  are  spaced  closely,  typically  no  greater  than  A/2  [8]. 
Accurate  performance  modeling  of  such  an  array  requires  proper  accounting  of  the  acoustic 
interaction  between  array  elements. 

Historically,  the  Pritchard’s  approximation  [9]  has  been  employed  to  evaluate  inter-element  ar¬ 
ray  interactions.  It  has  been  employed  to  model  a  bi-planar  array  by  the  Naval  Underwater 
Systems  Center  (NUSC,  New  London,  CT;  now  the  Naval  Undersea  Warfare  Center,  NUWC, 
Newport,  RI)  to  predict  the  accurate  beam  pattern.  Beam  pattern  predictions  of  computer  mod¬ 
els  were  then  compared  against  results  of  deep  water  testing.  Although  the  computer  models 
resulted  in  much  better  beam  patterns  with  regard  to  the  non-interacting  point  sources  model,  it 
showed  a  significant  deviation  at  and  near  end  fire  steering  angles  [8]. 
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Another  example  of  evaluating  array  interaction  comes  from  a  fundamental  acoustic  course  at 
the  Naval  Postgraduate  School  (NPS).  It  is  a  three-element  array  with  omnidirectional  drivers 
that  was  investigated  for  several  different  configurations  in  an  anechoic  chamber.  The  investiga¬ 
tion  uncovered  substantial  deviation  between  theoretical  and  experimental  results.  In  particular, 
near  end  fire,  source-level  amplitude  data  does  not  match  theory,  and  the  product  theorem  [10] 
cannot  explain  this  discrepancy  Figure  1.2. 


Polar  Plot  for 
847  Hz 


EXP. 

THEO. 
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Figure  1.2:  The  far-field  beam  pattern  for  a  three-element  array  with  omnidirectional  drivers  [3] 

In  the  early  1990s,  a  "T-Matrix  Method"  was  developed  by  NPS  researchers  Prof.  S.  R.  Baker 
and  Prof.  C.  L.  Scandrett  [6, 1 1]  in  order  to  properly  model  randomly  and  densely  packed  arrays 
of  active  sonar  transducers. 

To  build  an  array  model  using  the  T-Matrix  method,  one  must  obtain  the  free-field  radiation 
properties  for  a  given  electrical  drive  condition  and  the  scattering  properties  (T-Matrix  coeffi¬ 
cients)  for  a  single-transducer,  from  such  as  a  finite-element  model  (FEM). 

This  thesis  focuses  on  building  a  three-dimensional  (3D)  FEM  of  a  volumetric  piezoelectric 
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transducer  for  T-matrix  calculation  in  COMSOL  and  properly  developing  a  model  of  a  closely 
packed  linear  array  by  using  the  T-Matrix  Method. 

1.1.1  Related  Works 

Although  the  T-matrix  method  was  first  developed  in  the  1990s,  it  has  not  been  applied  to  a 
realistic  volumetric  array  model  because,  until  now,  the  capabilities  and  user-friendliness  of  the 
finite-element  computer  codes  available  for  performing  the  single-transducer  radiation  and  scat¬ 
tering  computations  required  for  the  T-matrix  method  were  insufficient.  Thus,  the  method  was 
not  of  practical  use  to  a  non-specialist  engineer.  The  attempt  of  prior  theses  to  compute  T-matrix 
elements  using  ATILA  finite-element  code  resulted  in  excessive  discrepancies  in  comparison  to 
analytical  results  [12, 13]. 

Recently,  that  situation  has  changed.  Using  the  COMSOL  multiphysics  finite-element  computer 
code,  finite-element  model  computations  for  the  acoustic  scattering  T-matrix  coefficients  of  a 
thin-wall  elastic  spherical  shell  now  show  excellent  agreement  with  those  computed  using  full 
3D  elasticity  theory  [14]. 

Therefore,  the  T-matrix  method  is  applicable  and  can  be  used  to  build  an  exact  model  of  densely 
and  randomly-packed  active  sonar  transducer  array.  Furthermore,  this  method  can  achieve  a 
more  accurate  analysis  of  array  performance  that  better  defines  the  acoustic  characteristics  of 
the  array. 

1.2  Purpose 

The  purpose  of  the  research  described  in  this  thesis  is  to  demonstrate  for  the  first  time  the 
practical  application  of  the  T-matrix  method  to  the  numerical  performance  modeling  of  an  ac¬ 
tive  sonar  array,  that  is,  to  combine  the  results  of  finite-element  model  computations  of  the 
single-transducer  free-field  radiation  and  scattering  characteristics  with  an  analytic  model  for 
the  free-field  propagation  between  transducers.  Specifically,  the  T-matrix  method  will  be  ap¬ 
plied  to  model  the  source  level  and  far-field  beam  pattern  for  several  example  unsteered  linear 
active  sonar  arrays  comprised  of  either  1)  radially-polarized  piezoelectric  spherical  thin-shell 
transducers  or  2)  type  V  flextensional  "ring-shell"  transducers,  model  34A0610,  manufactured 
by  Sparton  of  Canada  [7, 15].  The  former  was  chosen  as  an  academic  example,  as  approx¬ 
imate  analytic  radiation  and  scattering  properties  for  such  transducers  are  known;  the  latter 
was  chosen  because  it  represents  a  real  high-power,  low-frequency  sonar  projector  which  was 
built  and  tested  in  a  24-element,  densely-packed,  bi-planar,  active  array,  and  for  which  single- 
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element  and  array  measured  performance  characteristics  are  available  [7, 16].  Although  only 
unsteered  linear  array  configurations  of  1/2- wavelength  spacing  and  "ka"~l  were  considered 
herein,  the  ultimate  goal  of  the  research  program,  of  which  the  present  effort  is  a  contribu¬ 
tion,  is  to  compare  the  characteristics  of  the  as-tested  bi-planar  array  as  computed  using  the 
T-matrix  method  to  both  its  measured  characteristics  and  its  characteristics  as  computed  using 
the  so-called  "Pritchard  approximation"  to  model  array  element  interaction  [9]. 

1.3  Organization 

The  remainder  of  the  thesis  consists  of  six  chapters. 

Chapter  2  clarifies  features  of  COMSOL’s  Acoustic  Module  and  presents  the  addition  theorem 
used  in  coordinate  translation  matrices  for  an  array.  The  second  part  of  chapter  2  describes  the 
far-field  beam  pattern  calculation  for  a  linear  array  using  the  T-matrix  method. 

Chapter  3  describes  the  characteristic  features  of  the  modeled  transducers  which  consist  of  a 
radially  polarized  spherical  transducer  and  a  type  V  Flextensional  Ring-Shell  Model  34A0610. 

Chapter  4  is  devoted  to  building  3D  models  of  the  radially  polarized  piezoelectric  spherical  thin- 
shell  transducer  and  the  type  V  flextensional  ring- shell  transducer.  The  second  part  of  Chapter 
4  includes  an  evaluation  of  sought-after  transducer  parameters  from  the  3D  FEM  model. 

Chapter  5  calculates  the  total  acoustic  pressure  field  and  produces  the  far-field  beam  pattern 
plots  for  equally  spaced  linear  arrays  using  the  T-matrix  methodology. 

Chapter  6  concludes  this  thesis  and  provides  suggestions  for  future  work. 

1.4  Assumptions 

Single-frequency  harmonic  time  dependence  of  the  form  elC0t  is  assumed  throughout  this  thesis. 
The  symbol  is  used  for  the  imaginary  unit  to  avoid  confusion  with  the  spherical  Bessel 
function  of  the  first  kind,  jn. 
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CHAPTER  2: 
Theory 


2.1  COMSOL  Finite  Element  Analysis 

2.1.1  COMSOL  Multiphysics 

The  COMSOL  Multiphysics  numerical  package  is  a  finite  element  code  that  includes  simula¬ 
tion  software  capable  of  modeling  many  different  physical  phenomena  employing  a  variety  of 
different  physics-based  modules.  The  modules  involve  a  solution  method  of  Partial  Differential 
Equations  (PDE)  and  predefined  multiphysics  interfaces.  Predefined  multiphysics  interfaces 
simplify  the  model  by  using  a  built-in  coupled  physics  node.  They  undertake  the  analysis  for 
different  types  of  studies  (e.g.,  modal,  eigenfrequency,  frequency  response)  [17]. 


2.1.2  The  Acoustic-Piezoelectric  Interaction  Interface 

The  Acoustic-Piezoelectric  Interaction  Interface  is  one  of  the  physical  modes  of  COMSOL 
Multiphysics  within  the  Acoustics  Module.  This  module  can  solve  the  far-field  or  near-field  ra¬ 
diation  problem  by  incorporating  acoustic-structure  interactions  of  piezoelectric  materials  [18]. 
The  boundary  conditions  include  the  pressure  load  from  fluid  to  solid  and  a  structural  accelera¬ 
tion  between  the  solid  and  fluid.  The  Acoustic-Piezoelectric  Interaction  Interface  combines  the 
piezoelectric  device,  pressure  acoustics,  frequency  domain,  solid  mechanics,  and  electrostatic 
interfaces  [18], 


2.1.3  Acoustics  Interface  Theory  in  COMSOL  Multiphysics 

The  acoustic  wave  equation  is  expressed  in  COMSOL  by  the  following  equation  [18]: 

-^^+v'(Tv'’~q,'))=e”  <2'l> 

Bulk 

Modules 
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Variable 

Description 

Unit 

P 

Fluid  Density 

kg/m3 

c 

Speed  of  Sound 

m/s 

p 

Pressure 

Pa 

Qd 

Dipole  Source 

N/m3 

Qm 

Monopole  Source 

1/s2 

In  the  case  of  plane  waves,  the  pressure  varies  with  time  as 

P(x,t)  =p(x)el(0t.  (2.2) 

The  assumption  of  the  frequency  domain  or  time-harmonic  dependence  for  the  acoustic  wave 
equation  results  in  the  inhomogeneous  Helmholtz  equation  [18]: 

V'(-l(Vp-  <K))-^  =  a»  (23) 

Boundary  conditions  are  required  to  solve  the  inhomogeneous  Helmholtz  equation.  COMSOL 
provides  several  different  types  of  boundary  conditions,  depending  upon  the  analysis  performed. 
The  prescribed  boundaries  between  the  fluid-solid  boundary  are  added  after  nodes  are  selected 
and  fixed.  Three  built-in  boundary  conditions  are  available  [18]  on  the  boundaries  between  the 
fluid  and  solid: 

1.  free 

2.  sound  hard  boundary(wall) 

3.  acoustic  structure  boundary 

A  free  boundary  condition  assumes  the  medium  reacts  without  any  constraints  acting  on  the 
boundary. 

A  sound  hard  boundary  is  added  for  a  rigid  boundary.  With  this  boundary  condition,  the  normal 
component  of  differential  velocity  is  zero  [18].  This  boundary  condition  is  also  named  the 
Neumann  condition: 

n- (^(vP-qd))  =°  (2.4) 

An  acoustic  structure  boundary  couples  the  fluid  load  and  structural  acceleration  at  the  fluid-solid 


6 


boundary.  It  cannot  be  added  to  other  boundaries  [18]. 


F  P  =  -n p 
an  =  n  •  u„ 


(2.5) 


Variable  Description 

F  Force 

p  Pressure 

n  Unit  normal, outward  pointing  from  solid  domain 

an  Normal  acceleration  for  the  acoustic  fluid 

utt  Second  derivative  of  the  Structural  Displacement  with  respect  to  time 


2.1.4  Piezoelectric  Device  Theory  in  COMSOL  Multiphysics 

Piezoelectric  materials  are  solids  for  which  the  mechanical  strain  and  the  electric  field  in  the 
solid  influence  each  other.  Changes  in  strain  or  electric  field  cause  responses  of  the  polarization 
and  stress  within  the  solid.  For  example,  an  applied  potential  difference  induces  a  deformation 
of  the  piezoelectric  device,  and/or  the  deformation  of  the  piezoelectric  material  induces  an 
electrical  potential  difference  within  the  piezoelectric  device.  The  piezoelectric  relations  for 
linear  systems  having  small  strains  are  expressed  by  the  following  equations  in  COMSOL  [18]. 

Strain-Charge  : 

S  =  sET  +  dTE  (2.6) 

D  =  dT  +  £TE  (2.7) 

This  relation  can  be  rewritten  into  stress-charge  form  as 
Stress-Charge: 

T  =  cES  -  eTE  (2.8) 

D  =  eS  +  £5E  (2.9) 

Piezoelectric  ceramics  are  “transverse  isotropic”  materials  [19].  In  material  science,  orthogonal 
axes,  designated  “1,”  “2,”  and  “3”  are  defined  by  convention  for  a  given  material  type.  The 
polarization  is  conventionally  designated  as  the  “3”  direction,  and  COMSOL  uses  the  same 
convention  for  piezoelectric  devices. 
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The  piezoelectric  material  properties  are  expressed  in  stiffness,  compliance,  coupling,  and  di¬ 
electric  matrices  with  respect  to  their  local  coordinate  system.  The  default  local  coordinate 
system  corresponds  to  the  Cartesian  coordinate  system  [17]. 


Variable 

Description 

T 

Stress 

S 

Strain 

D 

Electrical  displacement 

Elastic  compliance  constant 

ce 

Elastic  stiffness  constant 

£s 

Dielectric  constant 

e 

Piezoelectric  stress  constant 

d 

Piezoelectric  strain  constant 

2.2  Piezoelectric  Spherical  Thin-Shell  Transducer 

The  spherical  thin  shell  transducer  is  a  particular  transducer  in  which  the  spherically  symmetric 
mode  (lumped  constant  mode)  can  be  developed  from  an  equation  of  motion  rather  than  a  wave 
equation  (Figure  2.1)  [4,20],  Therefore,  fundamental  acoustic  features  (e.g,  mechanical  Q, 
radiation  load)  of  a  thin  shell  transducer  can  be  accurately  represented  by  a  simple  closed  form 
expression.  The  equivalent  circuit  [4]  represents  a  model  of  a  thin  shell  transducer  (Figure  2.2). 
Hence,  it  is  particularly  convenient  for  illustrating  equivalent  circuits  and  design  principles  [20], 

A  radially  polarized  spherical  thin  shell  transducer  of  mean  radius  “a”,  and  thickness  “f”  is 
modeled.  The  radial  direction  coincides  with  the  material  ”3“  direction.  Electrodes  are  located 
on  the  inner  and  outer  surfaces  of  the  sphere  to  form  equipotential  surfaces.  Hence,  transverse 
fields  E\  and  Ei,  which  are  perpendicular  to  the  radial  direction,  are  neglected.  All  of  the  tensile 
stresses  are  equal  to  zero  except  the  circumferential  tensile  stress  T\  [4, 19].  The  relationship 
between  strain  and  stress  forces  in  the  “1”  and  “2”  directions  are  obtained  from  symmetry 
The  relationship  between  strain  and  stress  forces  in  the  1  and  2  directions  are  obtained  from 
symmetry 

T\  =  T2  =  tc 
S  i  =So=  sc. 

Compliance  coefficients  are  combined  to  form  a  circumferential  compliance  coefficient  ex- 
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O  V  Q 


Figure  2.1:  Radial  cross  section  of  a  radially  polarized  piezoelectric  spherical  thin  shell 
transducer 


pressed  as 

Sc  =  2  +512)- 

(2.10) 

The  piezoelectric  equations  are  then 

Sc  =  2s^Tc  +  d3\E3 

(2.11) 

D3  =  2d\  1  Tc  +  £33^3 . 

(2.12) 

The  electrical  displacement  equation  can  be  expressed  by  substituting  in  equation  (2.1 1) 


D3  = 


“1 r^c  +  e33  ( 1  —  kp2)  ^3  • 


(2.13) 


Where 


,  2_^31 


£33Sc 


(2.14) 


The  circumferential  strain  equals  the  ratio  of  the  radial  displacement  to  mean  radius  of  the 
sphere 

S,  =  -•  (2.15) 

a 


The  radial  displacement  can  be  obtained  from  the  equation  of  motion,  which  is  derived  from 
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Newton’s  Second  Law 


where, 


The  equation  of  motion  becomes 


C 

Fr  +  pwtaS  6  =  0 

otz 

Fr  =  F50  =  2Tcwt8Q. 


c  d2t 

pwtadG- =  —2  Tcwtdd. 
dtz 


(2.16) 


The  equation  of  motion  can  be  expressed  using  differentials  by  substituting  Equations  (2.11) 
and  (2.15)  into  Equation  (2.16) 


<92<§  1  e  J31 


aPsc 


dt 2  a2pSc 

The  free  resonance  frequency  of  the  undriven  sphere  (E=0)  can  be  written  from  2.17 

2  _  1 

=  a^' 

The  spherical  shell  displacement  for  harmonic  excitation  is  then  expressed  as 

■«2z7  „J(0t  „J_.,r.v2i 


(2.17) 


(2.18) 


ad^\  co^Eoe1®1  ad^iCO^Ei, 


CQq  —  coz 


(O^-CO2 


(2.19) 


After  using  the  relationship  between  symmetrical  strain  and  spherical  shell  displacement,  Sc  =  £,/ a, 
this  equation  can  be  written  in  a  more  convenient  form  for  computing  the  current,  I  =  4 na2£> , 


D3=£3(tI+^(i-^) 


SO 


I  —  4-7ta2i(oE3  +  e33  (l  —  ^p2)  )  • 

The  electrical  admittance  can  be  obtained  from  equation  (2.21) 

I  1  4 7ta2i(0  ( J31 


(2.20) 


(2.21) 


(On 


V  Z 


(ft)2-©2) 


+4(i-^2) 


(2.22) 
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The  equivalent  circuit  for  a  spherical  ring,  Figure  2.2,  can  be  obtained  for  the  electrical  admit¬ 
tance. 


Figure  2.2:  Equivalent  circuit  for  a  spherical  transducer  [4] 


for  which  the  electrical  element  values  are: 


_  A-KCld^i 

M  =  4  na2tp 

=  We3r3(l— 


rE  =  ^L 
m  4m 


7  —  n  rS  ika 

Z,r  —  p0LZ  l+ika 

S  =  4  na2 


(2.23) 


Element 

Description 

G0 

Electrical  loss  conductance 

Co 

Open-circuit  capacitance 

N 

Electromechanical  turns  ratio 

CE 

Short-circuit  compliance 

M 

Lumped  mass 

R-m 

Mechanical  Resistance 

Zr 

Free-field  radiation  impedance 

2.3  Type  V  Flextensional  Ring-Shell  Transducer 

Ring-shell  flextensional  transducers  are  desirable  because  of  their  ability  to  operate  at  low  fre¬ 
quencies  while  generating  high  acoustic  power.  They  also  have  higher  operating  depth  than 
uncompensated  flextensionals  of  similar  size  and  weight  [2],  Furthermore,  they  provide  wide 
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bandwidth  and  low  mechanical  quality  factor,  Qm,  which  are  very  important  characteristics  for 
underwater  applications. 

A  composite  piezoelectric  ceramic  and  steel  ring  is  used  as  a  motor  element.  The  ring  is  sand¬ 
wiched  between  two  convex  thin  spherical  shell  caps  Figure  2.3.  The  small  motion  of  the  ring 
causes  a  relatively  large  volume  displacement  of  the  shell.  The  ring-shell  transducer  is  depth 
compensated  and  operates  at  low  frequencies. 

The  desired  resonance  frequency  can  be  obtained  by  adjusting  the  piezo-ring  diameter.  In  order 
to  compensate  for  the  depth  pressure,  a  water  bladder  is  used.  The  bladder  is  opened  to  the 
ocean  and  adjusts  the  internal  pressure  according  to  the  depth  pressure.  [15] 


Figure  2.3:  Cross  Sectional  View  of  Type  V  Flextensional  Ring-Shell  Transducer 


2.4  T-Matrix  Method 

2.4.1  Introduction 

Application  of  the  T-Matrix  method  requires  computation  of  the  single-transducer  scattering 
properties  for  an  arbitrary  incident  wave.  Since  the  incident  pressure  distribution  is  not  known  a 
priori,  it,  along  with  the  single-transducer  free-held  radiation,  is  parametrized  as  a  series  expan¬ 
sion  with  to-be-determined  coefficients  using  the  free-held  eigenfunctions  of  the  linear  wave 
equation  in  spherical  coordinates.  Outgoing  waves  (free-held  radiated  and  scattered)  are  rep¬ 
resented  using  outgoing  traveling  wave  solutions  (spherical  Hankel  functions),  while  incident 
waves  are  represented  using  standing  wave  solutions  (spherical  Bessel  functions).  Unnormal¬ 
ized  spherical  harmonic  functions  are  used  in  all  cases  as  an  angular  basis  set,  as  shown  in 
equations  (2.24)  and  (2.25) 

oo  n 

p‘(r,e,<l>)=  E  E  Kmjn(kr)a:(e,(P)  (2.24) 

n  m=—n 
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(2.25) 


00  n 

p’(r,e,4,)  =  £  £  PZjV(kr)Q’;(e,, p), 

n  m=—n 

where 


£i“(e,0)  = 

P™(  cosd)eimt 

Variable 

Description 

M,*) 

The  spherical  coordinates 

pi 

1  nm 

Incident  pressure  coefficient 

^nm 

Scattered  coefficient  from  the  obstacle 

jn(kr) 

nth  order  spherical  Bessel  function 

hf\kr) 

nth  order  spherical  Hankel  function  of  the  second  kind 

Unnormalized  spherical  harmonic  function 

P„m(  COS0) 

Associated  Legendre  function  of  the  first  kind 

To  simplify  further  development,  the  following  matrix  notation  will  be  employed.  The  spherical 
harmonic  component  coefficients  for  any  quantity  belonging  to  a  single  radiator  are  arranged  in 
a  column  matrix,  denoted  using  bold  capital  letter,  e.g.,  by  “flattening  out”  the  double  sequence 
over  n  and  m,  as  follows 


/  ps  \ 

M)0  1 

(  ps  ^ 

Ps 

M  — 1 

n 

ps 

r10 

= 

ps 

?h 

n 

v;  J 

v-  / 

The  one  dimensional  running  index  “z”  of  the  column  matrix  equals 


r\ 

1  +n  +n  +  m 


In  principle,  the  length  of  such  matrices  may  be  infinite;  however,  in  a  practical  computation, 
the  number  of  spherical  harmonics  kept  is  truncated  to  a  maximum  value  up  to  a  certain  degree, 
typically  the  second  degree  (quadrupolar,  n=2)  [6].  If  the  maximum  spherical  harmonic  degree 
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2 

kept  is  nmax,  the  total  number  of  spherical  harmonic  components  kept  is  N  =  (nmax  +  1)  ,  e.g., 
9  for  nmax  =  2. 


2.4.2  Unnormalized  Versus  Normalized  Spherical  Harmonics 

The  unnormalized  spherical  harmonic  function  Cl™  (0,0)  is  related  to  the  normalized  spherical 
harmonic  function  Y™  (0, 0)  [21]; 

Y™  (0,0)=  CnmCl™  (0,0)  =  CnmP™  (cos  0)  exp  (/m0)  (2.26) 


where 


C„ 


(-ir 


2n  +  1  (n  —  or) ! 
47T  («  +  /«)! 


(2.27) 


Any  function  can  be  decomposed  into  a  superposition  of  either  normalized  or  unnormalized 
spherical  harmonic  functions  [21] 

oo  n 

/(0,0)=£  I  anmY™(d^)  (2.28) 

n=0m=—n 

where  the  constant  anm  is  given  by 

jr  2/r 

a,™  =  j  J  [y*  (  0 , 0 ) ]  V  ( 0 , 0 )  ^  =  C«m4m  (2.29) 

0  0 


so 

oo  ft 

/(»,*)=  E  L  (2.30) 

ft=0 m=—n 

The  orthonormality  property  of  the  associated  Legendre  equation  is  used  for  the  limit  0  to  n 
0-dependent  and  0  to  2 n  0-dependent. 

The  Inm  is  given  by  [14] 

71  271 

Inm  =  j  I  [n-(0,0)]7(0,0)JQ.  (2.31) 

0  0 
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2.4.3  T-Matrix  Calculation 


The  multiple  scattering  of  acoustic  waves  by  transducers  in  an  array,  to  all  orders,  is  rigor¬ 
ously  and  self-consistently  accounted  for  by  enforcing  that  the  scattered  component  amplitudes 
at  each  transducer  are  related  to  the  free-field  incident  component  amplitudes  from  all  other 
transducers  by  a  transition,  or  T-,  matrix 


ps  =  [j]  pi 


(2.32) 


Employing  the  matrix  representation  introduced  in  Section  2.4.1,  this  relationship  can  be  ex¬ 
pressed  as 


'  Pf ' 

'  Tn  . 

'  P\  ' 

Pfr 

.  Tnn 

Ph 

(2.33) 


2 

where  N  =  (nmax  +  1)  • 

The  T-matrix  elements  are  computed  through  a  succession  of  scattering  computations,  as  fol¬ 
lows.  An  incident  wave  is  imposed  with  unit  amplitude  in  exactly  one  member  (  component)  of 
the  incident  eigenfunction  expansion,  and  the  amplitudes  of  all  members  (components)  of  the 
scattered  eigenfunction  expansion  found.  This  determines  the  member  of  exactly  one  column 
of  the  T-matrix.  The  T-matrix  elements  of  that  column  equal  the  component  pressure  ampli¬ 
tudes.  For  example,  to  determine  the  first  column  of  the  T-matrix,  elements  denoted  Rnm  below, 
an  incident  wave  with  /)')(J  =  1  and  all  other  P'nm  =  0  is  imposed,  and  the  resulting  left  hand 
side  column  matrix  of  scattered  wave  component  coefficients  P%m  determined  by  applying  or¬ 
thogonality  of  the  unnormalized  spherical  harmonic  functions  (see  Section  2.4.2  for  details).  If 
ps  (a,  6,  (f>)  is  the  scattered  pressure  amplitude  at  radius  a  in  the  fluid,  then 


where 


C2  I 

h (, ka ) 
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inm  =  j  I  [n™(0,0)]V(fl,0,^)dQ. 

0  0 


(2.34) 


(2.35) 
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ps 

r00 


ps 

1  nm 


Roo 


Rn 


where 


0 


(2.36) 


Rnm  —  Pnm  — 


C2  I 

_ nm 1n 

h $  ( ka ) 


(2.37) 


This  follows  from  setting 


oo  n  oo  n 

E  E  O>i2)(fa>«(0,<P)  =  E  E  ct,/™Q7(e,*) 

w  m=—n  n=0m=—n 


(2.38) 


It  should  be  noted  that  by  its  spherical  symmetry  property,  the  T-matrix  for  the  spherical  thin- 
shell  transducer  must  be  diagonal.  Hence,  the  scattering  cannot  mix  components.  This  is  not 
exactly  the  case  for  the  ring-shell  transducer.  By  its  axial  symmetry  property,  scattering  by  the 
ring-shell  transducer  can  mix  “n”  values,  but  not  “m”  values. 


2.4.4  The  Radiation  Pressure  Coefficient  Calculation 

In  the  absence  of  scattering,  the  free-held  radiated  pressure  is  expressed  with  the  spherical 
Hankel  function  and  unnormalized  spherical  harmonic  function 

oo  n 

P*(r,0,t)='E  E  rfmh*'  (kr)  Cl™  (0,  *)  .  (2.39) 

n=Qm=—n 


The  radiated  pressure  coefficients  can  be  computed  using  the  FEM.  After  comparing  the  spher¬ 
ical  harmonic  function  expansion  (2.28),  the  radiation  pressure  amplitude  can  be  expressed  in 
terms  of  Inm  as 


C2  1 

pR  _  ^ nm1nm 


h ^  ( ka ) 


(2.40) 


The  right-hand  side  variable  seems  the  same  as  the  scattered  pressure  amplitude,  but  the  surface 
integral  (Jnm)  in  this  case  will  be  different  as  scattering  has  been  neglected. 
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2.4.5  The  Addition  Theorem 


In  the  case  of  an  array,  a  single  common  coordinate  system  is  inconvenient  in  applying  boundary 
conditions.  The  spherical  addition  theorem  must  be  used  in  this  situation.  It  provides  for  a 
representation  of  a  pressure  field  relative  to  one  coordinate  system  in  terms  of  another. 

Consider  two  active  spherical  sources  of  radius  “a”  and  separated  by  distance  “d”,  as  shown 
in  Figure  2.4.  The  outgoing  pressure  of  sphere  two  in  mode  ( n ,  m )  can  be  determined  as  a 
superposition  of  free-field  incident  pressures  on  sphere  one  in  modes  (s,  u): 

O®  (kr2 )  C  (ft ,  fc )  =  E  E fij': ’  (*»)  • K  (e. .  *l )  ,  <2.41 ) 

5  u 

where 

(K,u)  =  [G12]  (O  (2-42) 

and  G\2  can  be  written  from  [5] 

n-\-s 

E  a(^(,n,».m)/.!2,(«rr12)£ir“(ei2.'fl2)  .  (2.43) 

t  =  \n  —  s\ 
t  >  \m  —  u\ 


where 


Y1 


Figure  2.4:  Two  coordinate  frames  with  arbitrary  relative  positions  and  orientations  from  [5] 
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a(s,t,n,u,m)  =  is+t  n 


(2s  +  l)  (2t  +  l)  (s  —  u)\  (t  —  m  +  u)  \  (n  +  m)\[(s  +  t  +  n)  /2]  \  1 
[(n-H-s)  /2}\  [(n  +  s-t)  /2]\  [(s  + 1  -  n)  /2]l  (s  + 1  +  n  +  l)  \  J 


w  max 

E  (-i> 

w=wmin 


s  +  t  —  n 

(s  +  t  —  n)  /2  +  w 


(S  +  tl  —  t 

T  n  —  r )  j2  —  u  +  xv 


t  H  ~\ - S 

(t  +  n  —  s)  /2  +  m  —  u  +  w 


with 


a 

b 


a\ 

&!(«-&)! 


2wmin  =  max  {n  —  s  —  t,s  —  n  —  t  —  2m  +  2  u,t  —  s  —  n  +  2  u) 


2wmax  =  min  (s  +  s  —  n,n  +  t  —  s  —  2m  +  2 u,  n  +  s  —  t  +  2 u)  . 


Variable  Description 

r\2  the  inter-element  spacing  distance, “rf.” 

012  the  polar  angle  of  origin  “1”  with  respect  to  origin  “2” 

012  the  azimuthal  angle  of  origin  “1”  with  respect  to  origin  “2” 

For  a  line  array  along  the  polar  axis,  which  is  all  that  will  be  considered  herein,  0i2  equals  either 
0  or  71,  and  0 12  is  “0”. 

2.4.6  Total  Acoustic  Pressure  Calculation 

The  multi-scattering  of  acoustic  waves  must  be  considered  in  the  acoustic  field  calculation.  The 
acoustic  field  is  regarded  as  a  superposition  of  outgoing  waves  from  each  transducer  in  the 
sonar  array  using  the  eigenfunctions  of  the  wave  equation  in  spherical  coordinates  as  a  basis  set 
(e.g.,  spherical  Hankel  functions  and  spherical  harmonics)  centered  at  each  array  element.  Out¬ 
going  wave  components  account  for  both  the  free  field  radiation,  when  a  transducer  is  electri¬ 
cally  energized,  and  the  scattered  field. 
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pout  _  PR  +  PS 


(2.44) 


“Self  consistency”  is  imposed,  i.e.,  the  incident  pressure  on  one  transducer  is  due  to  the  outgoing 
pressure  from  all  other  transducers.  This  leads  to  a  coupled  set  of  linear  equations  in  terms  of 
the  radiated  and  scattered  pressure  coefficients.  For  example,  for  two  elements 


(p?)  =[T1Gu]  (l?  +  Pf) 

and 

(pf)  =[t2g21]  (pf+pf)  . 

Variable  Description 

Scattering  pressure  amplitude  for  Ith  transducer 
G iq  Translation  formula  from  system  q  to  system  / 

T;  Scattering  pressure  amplitude  for  Ith  transducer 
Pf  Free-field  radiation  pressure  amplitude  for  Ith  transducer 


In  the  matrix  form  (2.45)  and  (2.46)  can  be  expressed  as 


I  TrG  12 

T2G2i  I 


pS 

M  1  - 

pS 
r2 


TiG12P« 

T2G21Pf 


a  = 


The  “a”  matrix  is  defined  as 

0  T!G12 

T2G21  0 

After  substitution  into  (2.47)  we  get  the  simpler  matrix  form: 


[I-al 


pS 

r2 


1  '  =  a 


Pf 

Pf 


or  equivalently 


pS  )  =[!-«]  '« 

r2 


pf 

pf 


This  result  is  easily  generalized  to  the  case  of  “N”  transducers.  The  result  is 


(2.45) 

(2.46) 


(2.47) 


(2.48) 


(2.49) 


(2.50) 
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(  Pls  \ 


/ 


(I -a)  !a 


/  P?  \ 


VPN  / 


From  this  equation  the  scattering  pressure  amplitudes  can  be  found. 

The  summation  of  the  scattered  and  radiated  pressure  for  each  individual  coordinate  system  of 
the  transducer  provides  the  total  pressure  at  any  field  point  with  respect  to  a  global  origin  in 
spherical  coordinates. 


A^X  %nax  n 

pM,*)  =  E £  £  (2-51) 

/=]  n= 0  m=-n 

where  r/,  0/,  and  (f>/  are  the  range,  polar  angle,  and  azimuth  angle,  respectively,  of  the  field  point 
with  respect  to  the  local  coordinate  axes  of  transducer  In  all  of  the  arrays  considered  herein, 
the  local  coordinate  axes  of  all  the  transducers  are  assumed  the  same;  only  the  origins  differ  in 
location.  Hence,  0/  =  0  and  <j>i  =  <j>  for  all 


2.4.7  Far  Field  Beam  Pattern 

A  sonar  array  can  be  represented  as  a  collection  of  arbitrarily-randomly  and  densely-packed 
acoustic  transducers.  An  example  for  a  sonar  array  is  shown  in  Figure  2.5,  a  bi-linear  array 
spherical  thin  shell  transducer  elements.  The  horizontal  far-field  beam  pattern  of  this  sonar  array 
is  represented  in  Figure  2.6  [6],  The  amplitude  of  pressure  is  given  in  terms  of  a  directional 
factor,  and  the  pressure  amplitude  along  the  acoustic  axis  of  source,  Pax, 

P(r,e,<t>)  =  Pax(r)H(0,(t>)  .  (2.52) 


The  pressure  amplitude  change  is  proportional  to  1  /r  in  the  far  field,  so  equation  (2.52)  becomes 


lim  p  (r,  0,0) 

r— oo 


Pax(l)H(e,  0) 
r 


(2.53) 


where 

Pax{  1)  =  lim  p  (r,0,0)  r 

r— o 

is  the  far-field,  on-axis  pressure  amplitude,  “referred”  to  a  range  of  1  m  . 
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ORIENTATION  OF  SPHERICAL  TRANSDUCERS 


Figure  2.5:  Orientation  of  the  spherical  transducers  in  the  16  element  array 


Equation  (2.51)  can  be  simplified  using  the  asymptotic  form  of  the  spherical  Hankel  function, 
which  can  be  derived  from  the  standard  Hankel  function  [22] 

h")(kr)=\iwrHnk(kr)-  {2M) 

In  the  far  field,  the  asymptotic  form  of  the  Hankel  function  is  inversely  proportional  to  the 
square  root  of  the  distance 


lim  if  ,  (kr)  =  w 

r^-oo  n+2  V  71  kr 


(2)  (b-A  _  .  /  ^  „-ikr  -i(-§ («+2)-f) 


(2.55) 


So,  the  asymptotic  spherical  Hankel  function  can  be  written  as 


r— o 


71 


lim  hP  (kr)  =  J—*I  —  ) 


2 kr  V  7tkr 


After  simplification  it  becomes 


lim  h ^  (kr)  =  — 

r— AT 


(2.56) 
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Figure  2.6:  Source  level  of  array,  ka=l,  analytically  scaled  by  far-field  of  a  single  spherical 
transducer  [6] 


Applying  this  form  to  equation  (2.51),  and  noting  that  “kr”  above  represents  k  •  r,  the  total 
acoustic  pressure  for  the  far-field  limit  can  be  obtained 


NX  e—ik-(r—?0l)  "max  "  r 

Pff(r,0,<l>)  =  limp(r,0,<j>)  =  '£r  z  £  i(w+1)  £  Pimn  +  rfmn  , 

r^°°  /=1  k-yr  —  rgi)  n=0  m=-nl  J 

(2.57) 

where  the  direction  of  k  equals  the  direction  of  r  and  r(>/  is  the  vector  position  of  the  local  origin 
of  axes  for  transducer  with  respect  to  the  global  origin. 


In  the  far  field,  all  radial  vectors  ( r  —  roi )  are  approximately  parallel.  The  distance  to  the  far 
field  does  not  vary  much  for  each  element,  so  the  factor  roi  can  be  neglected  compared  to  r 
in  the  amplitude  factors  (denominators).  Although  the  path  length  difference  can  be  ignored 

for  amplitude,  it  has  to  be  taken  into  account  for  phase.  Factoring  out  the  common  phase  —  kr, 

— > 

the  factor  k  ■  roi  is  the  phase  difference  between  the  global  origin  and  the  acoustic  origin  of 
transducer  /.  If  the  acoustic  center  of  transducer  /  is  located  at  (. xoi,yoi,zoi), 


k-r0i=  kxx0i  +  kyy0i  +  k-z.oi 
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which  can  be  written  as 


k-foi  =  k (. xoi  sin  6 cos </>  +yoi  sin  Q  sin <j>  +  zoi cos  0)  •  (2.58) 

After  the  substitution  of  this  equation  into  equation  (2.57),  the  far-field  total  acoustic  pressure 
equation  is  obtained 


y—ikr  Nx  _ 


kr 


Jk-rgi 


1=1 


klmsix 

E< 

n= 0 


(»+!) 


^Imn  +  ^Jn”(9,«.)  (2.59) 


with  k  •  7qi  given  by  (2.58) . 
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CHAPTER  3: 

Spherical  and  Ring-Shell  Transducer  Descriptions 


The  radially  polarized  spherical  thin-shell  transducer  and  the  class  V  flextensional  ring-shell 
transducer  are  picked  as  example  applications  for  the  T-matrix  calculation  of  a  volumetric  array 
transducer. 


3.1  Piezoelectric  Spherical  Thin-Shell  Transducer 

The  piezoelectric  ceramic  thin  shell  is  polarized  in  thickness  (radial)  direction,  and  electrodes 
are  installed  on  the  inner  and  outer  surface  of  the  thin  shell  to  provide  radial  polarization.  The 
shell  is  made  of  lead  zirconate  titanate  (PZT-8).  This  material  is  treated  as  a  member  of  the 
Navy  Type  III.  The  Navy  Type  III  is  defined  as  “hard,”  which  delineates  a  relative  difficulty 
of  depolarization  of  the  ceramic  [19].  PZT-8  is  preferred  for  high  powered  transducers.  The 
properties  of  PZT-8  are  obtained  from  the  COMSOL  material  library. 

The  thin  shell  is  surrounded  by  a  water.  The  inner  domain  of  the  thin  shell  is  defined  as  air. 
The  frequency  of  the  transducer  is  defined  with  respect  to  “ka=l”  which  is  the  typical  “ka” 
value  of  low  frequency  transducers  [23].  Main  operational  and  dimensional  parameters  are  the 
following: 


Parameter 

Value 

Notation 

Frequency 

474  Hz 

/ 

Applied  Electric  Potential 

1  V 

V 

Shell  outer  Radius 

0.5  m 

a 

Shell  thickness 

0.01  m 

t 

Radiation  boundary  dimension 

2.5  m 

R 

Table  3.1:  The  spherical  thin-shell  transducer  dimensional  and  operational  parameters 


The  thin-shell  transducer  can  be  properly  represented  by  an  equivalent  circuit  and  the  radiation 
load  can  be  calculated  by  again  using  an  equivalent  circuit  (See  for  details  Chapter  2.2). 


3.2  The  Class  V  Flextensional  Ring-Shell  Transducer 

The  class  V  flextensional  ring-shell  transducer  (Sparton  Model  34A0610),  is  shown  in  Fig¬ 
ure  3.1,  and  has  been  selected  as  a  realistic  model  of  a  volumetric  transducer.  The  ring-shell 
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transducer  is  versatile  in  the  case  of  low-frequency  and  high-power  applications  [2],  It  is  com¬ 
prised  of  a  composite  piezoelectric  ceramic  and  steel  ring  sandwiched  between  two  thin  convex¬ 
shaped  steel  shell  segments.  The  radial  motion  of  the  piezoelectric  ring  causes  large  volume 
displacement  over  the  shell  and  provides  high  acoustic  output  [15]. 


Figure  3.1:  The  depth-compensated  ring  shell  projector  [2] 

3.2.1  Characteristics  of  the  Transducer 

The  outer  piezo-ring  diameter  is  81  cm.  The  ring  is  made  of  piezoelectric  material  (PZT-4), 
which  is  a  motor  element  of  the  transducer.  PZT-4  is  classified  as  a  “hard”  ceramic  type  (Navy 
Type  I)  and  is  used  primarily  for  high-powered  transducers  [19].  The  piezo-ring  consists  of  144 
plates  of  thickness-poled  lead  zirconate  titanate,  ceramic  staves  connected  electrically  in  paral¬ 
lel.  Rectangular-shaped  ceramic  staves  are  separated  by  72  steel  wedges  to  create  an  annular 
shape  form  [7],  Each  steel  wedge  separates  two  bonded  ceramic  plates.  Steel  wedges  are  axially 
longer  than  ceramic  staves  in  order  to  protect  a  ceramic  surface  from  the  stress  that  is  caused  by 
the  shell  edge  [2].  A  tension  applied  (25-40  MPa)  fiber  wrapping  encloses  the  piezo-ring.  The 
fiberglass  wrapping  is  bonded  with  an  epoxy  resin.  The  fundamental  operational  characteristics 
of  transducer  are  given  by  [1,15, 24] 

The  piezoelectric  elements  are  polarized  tangentially  and  are  connected  in  parallel.  The  polar¬ 
ization  direction  is  shown  in  Figure  3.2. 

In  common  with  the  flextensional  transducer,  the  operating  depth  of  the  ring-shell  projector  is 
limited  to  a  few  meters  without  a  depth  compensating  system.  Different  systems  (e.g.,  liquid, 
compliant  tubes,  self-supporting,  etc.)  may  be  used  to  overcome  the  depth  limited  operation. 
The  class  V  flextensional  ring-shell  projector  is  depth  compensated  by  an  internal  water  bladder 
that  is  opened  to  the  sea  through  the  water  inlet  [15].  The  water  bladder  is  installed  into  the 
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Parameter 

Frequency 

Maximum  Source  Pressure  Level 

Efficiency 

Operational  Depth 

Bandwidth  at  Resonance(-3  dB  points) 


Value 

610  Hz 

213  dB  re  1  pPa  at  1  m  (driven  by  3000  Vrms) 

%90  at  resonance,  %65  at  -6  dB 

Exceeds  400  m 

>140  Hz  (470  Hz  to  750  Hz) 


Table  3.2:  Ring-Shell  Transducer  Operation  Characteristic 


interior  cavity  of  the  shell  segments  (Figure  3.3).  The  air  domain  between  the  shell  segments  is 
compressed  by  the  water  bladder  with  respect  to  the  depth  pressure.  The  water  bladder  depth 
compensation  system  decreases  the  damping  and  supports  shells  against  the  depth  pressure 
without  constraining  motion  of  the  shell. 


The  in-air  weight  without  a  foundation  is  245  kg,  while  the  in-water  weight  is  272  kg. 


Figure  3.3:  Cross  sectional  view  of  ring  shell  projector  [2] 
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CHAPTER  4: 
The  3D  Models 


The  numerical  finite-element  models  of  a  radially  polarized  thin-shell  spherical  transducer  and 
a  class  V  flextensional  ring-shell  transducer  were  developed  using  the  COMSOL  Multiphysics 
version  4.3a  finite-element  code  [25],  The  acoustic-piezoelectric  interaction ,  frequency  domain 
interface  of  the  acoustic  module  was  selected  using  multi-physics  because  the  fluid- structure 
coupling  properties  are  enhanced  at  boundaries  (See  Chapter  2.1.3  for  details). 

The  development  of  models  is  comprised  of  the  following  steps: 

1.  Selection  of  the  COMSOL  module  and  physics  interface 

2.  Adding  the  required  type  of  analysis  (i.e, frequency  domain,  eigenfrequency,  frequency- 
domain  modal) 

3.  Creating  a  geometry  in  the  selected  space  dimension  (i.e,  3D,  2D,  2D  axisymmetric) 

4.  Setting  material  properties 

5.  Assigning  material  geometrically 

6.  Applying  the  physics  interaction  to  the  model 

7.  Meshing 

8.  Computation 

9.  Postprocessing 

4.1  General  Model  Settings 

Development  settings  are  the  same  for  both  models. 
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1.  The  built-in  acoustic-piezoelectric  interaction, frequency  domain  interface  (ACPZ)  was 
used.  This  interface  couples  the  field  boundary  conditions  between  fluid  and  the  piezo¬ 
electric  material.  The  normal  component  of  acceleration  over  the  piezoelectric  material’s 
surface  generates  sound  within  the  fluid.  Inversely,  the  piezoelectric  material  undergoes 
a  mechanical  distortion  due  to  applied  fluid  acoustic  pressure  over  its  surface  [18]. 

2.  The  free-field  incident  waves  are  referred  to  as  “the  background  pressure  field ”  in  COM- 
SOL  [18].  The  background  pressure  field  is  defined  for  each  spherical  harmonic  compo¬ 
nent  (n,  m).  In  each  computational  run,  exactly  one  incident  component  is  applied  within 
the  fluid  domain.  The  incident  pressure  amplitude  was  set  to  “1”  so  the  incident  standing 
wave  equation  (2.24)  becomes 

oo  n 

P‘(r, e, *)=£  £  j„(kr)0»( 9,*). 

n  m=—n 

3.  The  spherical  wave  radiation  condition  is  applied  on  the  outer  boundary  of  the  fluid  do¬ 
main  at  radius  2.5  m  to  reduce  the  reflection  on  the  boundary  and  to  truncate  the  numerical 
domain.  It  models  acoustic  propagation  without  reflection  from  the  boundary. 

4.  The  reference  pressure  for  sound  pressure  level  (SL)  was  selected  as  water  (Ip  Pa)  from 
the  SL  settings  window  in  COMSOL. 

5.  The  maximum  size  of  mesh  elements  was  limited  with  respect  to  the  wavelength  of  sound 
in  water.  To  achieve  a  reliable  solution,  the  wavelength  has  to  be  resolved  by  a  10-12  de¬ 
grees  of  freedom  (DOFs)  corresponding  to  six  elements  in  a  wavelength.  In  other  words, 
the  mesh  element  size  should  be  no  larger  than  A/5  [18].  This  criteria  was  applied  strictly 
to  each  mesh  sequence. 

6.  The  T-matrix  elements  were  calculated  for  monopole,  dipole,  and  quadrupole  spherical 
harmonic  components. 

7.  Two  studies  have  been  examined  in  each  model,  i.e,  study  1  for  calculation  of  the 
single-element  free-field  radiated  pressure  while  study  2  corresponds  to  determination 
of  the  T-matrix  elements,  each  using  the  post  processing  feature. 
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System  Requirements 

The  Apple  iMac  MC813LL/A  employed  for  COMSOL  has  the  following  technical  specifica¬ 
tions: 

a.  2.7  Ghz  Intel  Core  i5  processor 

b.  32  GB  1333  MHz  DDR3  memory 

4.2  The  Spherical  Thin-Shell  Transducer 

The  main  purpose  for  building  a  model  of  the  radially  polarized  piezoelectric  spherical  thin-shell 
transducer  is  to  calculate  the  elements  of  its  T-matrix  and  the  coefficients  of  its  free-field  radi¬ 
ated  pressure.  The  validation  of  the  COMSOL  Multiphysics  ACPZ  interface  was  evaluated 
by  comparing  the  obtained  values  of  the  free-field  radiated  pressure  field  at  the  surface  of  the 
spherical  thin-shell  transducer  with  an  equivalent  circuit  model  of  a  spherical  transducer. 

4.2.1  Geometry 

The  thickness  of  the  model  shell  is  0.01  m  and  has  a  radius  of  0.5  m.  The  piezoelectric  ceramic 
is  surrounded  by  water  (Figure  4.1).  The  cavity  of  the  thin-shell  is  air  filled. 


Figure  4.1:  Pictorial  representation  of  the  3D  spherical  thin-shell  transducer  and  fluid  domain. 

The  non-dimensional  parameter  ka=  1  was  selected  to  simulate  the  typical  transducer  operating 
condition  and  the  frequency  was  adjusted  to  enforce  ka=  1.  The  model  parameters  can  be  seen 
in  Table  3.1. 
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The  surfaces  of  the  spherical  thin-shell  were  covered  with  electrodes.  The  inner  surface  of 
thin-shell  was  grounded  while  the  outer  surface  was  driven  by  a  1  V  peak  amplitude  potential 
(Figure  4.2). 


▲  1 


Figure  4.2:  Electric  potential  of  the  thin-shell 


4.2.2  Material 

Three  different  material  types  were  assigned  to  the  geometry. 

Lead  Zirconate  Titanate(PZT-8) 

Properties  of  PZT-8  were  assigned  to  the  spherical  shell.  It  is  a  very  common  piezoelectric 
material  also  known  as  a  Navy  Type  III  that  is  often  used  for  power  transducers.  This  material 
provides  a  maximum  strength  against  high  tensile  forces.  Material  properties  are  already  built 
into  the  COMSOL  Material  Library,  and  are  shown  in  Table  4.1. 

The  default  polarization  direction  is  defined  along  the  “z”  axis  in  COMSOL  [18],  but  in  the 
present  case,  the  polarization  direction  is  radial.  In  order  to  obtain  this  polarization  a  user- 
defined  base  vector  system  was  added  to  the  “definition”  branch  . 


32 


Parameter 


Value  Unit 


Notation 


Density 

7600 

kg/m3 

Compliance  Matrix 

sn 

11.5 

10“12m2/N 

s\2 

-3.7 

10-12m2/N 

s13 

-4.8 

10_12m2/N 

S  33 

13.5 

10_12m2/N 

^44 

31.9 

10_12m2/N 

S66 

30.4 

10_12m2/N 

Coupling  Matrix 

d-}\ 

-97 

10"12C/N 

^33 

225 

10"12C/N 

^15 

330 

io~12c/n 

^15 

330 

io~12c/n 

Relative  Permittivity 

eh/£o 

1290 

1000 

£0 

8.85 

io~12f/m 

Table  4.1:  PZT-8  material  properties  (From  the  COMSOL  material  library) 


The  spherical  coordinate  system  was  defined  using  the  unit  basis  vector  (i,  j,  k)  [21] 


xl  =  ©0  =  cos  0  cos  <j>i  +  cos  0  sin  (j) j  —  sin  0k 

x2  =  =  -sin^)i  +  cos^j  (4.1) 

x3  =  er  =  sin  0  cos  (j) i  +  sin  0  sin  (j) j  +  cos  0k  . 

After  the  spherical  coordinate  system  was  added  in  the  definition  branch,  it  was  selected  as  a 
coordinate  system  from  the  selection  menu  in  the  material  properties  window. 

The  X3  unit  vector  of  the  spherical  coordinate  system  shows  the  polarization  direction,  and  in 
this  model  it  is  referred  to  as  the  radial  direction. 
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Water 

Water  was  used  to  define  the  fluid  domain  which  surrounds  the  transducer.  The  material  prop¬ 
erties  were  assigned  as  shown  in  Table  4.2. 


Parameter  Value  Notation 

Density  1000  kg/m3  p 

Speed  of  Sound  1490  m/s  c 

Table  4.2:  Water  material  properties 


Air 

Air  was  added  to  define  the  interior  cavity  of  the  thin-shell.  The  material  properties  are  built 
into  COMSOL. 


4.2.3  Meshing 

User-controlled  meshing  was  used  for  the  sequence  type  to  control  and  optimize  the  numerical 
mesh.  The  “free  mesh”  algorithm  is  preferred  over  the  geometry.  The  spherical  transducer 
domain  and  fluid  domain  used  a  finer  mesh  size  employing  tetrahedral  elements.  The  entire 
model  consists  of  421,598  elements  (Figure  4.3).  The  average  mesh  quality  is  0.73  where  the 
quality  measure  is  a  number  between  0  and  1  (Figure  4.4).  An  acceptable  mesh  quality  is 
typically  larger  than  0.1  according  to  the  COMSOL  user’s  manual  [17]. 

According  to  the  COMSOL  user’s  manual,  the  maximum  element  size  of  mesh  should  be 
smaller  than  A/5  (0.628  m  in  water  at  474  Hz)  in  an  acoustical  model  [18].  We  have  cho¬ 
sen  the  maximum  element  size  of  the  mesh  to  be  A/1 1.5  or  0.273  m  which  ensures  accuracy. 
The  statistical  mesh  quality  for  boundaries  is  represented  in  Table  4.3 


Boundaries 

Element  Type 

Amount 

Avg.  Element  Quality 

Radiation 

Triangular 

2,852 

0.98 

Piezo  material  outer 

Triangular 

19,416 

0.99 

Piezo  material  inner 

Triangular 

18,504 

0.99 

Table  4.3:  The  meshing  properties  of  the  3D  spherical  transducer  model  with  a  “finer  mesh” 
size  for  all  domains. 
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yix 


Mesh 


Figure  4.4:  Mesh  quality  of  the  model 
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4.2.4  The  Free-Field  Radiated  Acoustic  Pressure  Field 

The  free-held  radiated  acoustic  pressure  held  is  plotted  versus  the  on-axis  distance  from  the 
thin-shell  transducer  surface.  Figures  4.5  and  4.6  show  the  spatial  free-held  radiated  acoustic 
held.  Pressure  magnitude,  refer  to  one  meter  was  found  as 


P(l)  =  r  p  =  9.76  Pa.m  . 
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Free-Field  Radiated  Acoustic  Pressure  (Pa) 


Figure  4.5:  The  free-held  radiated  acoustic  pressure  held  (magnitude)  vs  distance  from  the 
center  of  the  spherical  thin-shell  transducer 


freq(l)=474 

TOTAL  ACOUSTIC  PRESSURE  FIELD  MAGNITUDE  (Pa) 

A  19.519 


T  5.74X10'5 

Figure  4.6:  The  multi-slice  plot  of  the  free-held  acoustic  radiated  held  (magnitude)  in  “Pa”  for 
a  modeled  spherical  thin-shell  transducer 


4.2.5  Approximate  Analytical  Radiation  Pressure  Calculation  for  a  Thin- 
Shell  Transducer 


The  free-held  radiated  acoustic  pressure  can  be  expressed  by  the  product  of  the  mechanical 
radiation  impedance,  Zr,  and  the  surface  velocity,  u,  integrated  over  the  surface  of  the  transducer 


P  = 


Zru 
4  na2  ' 


(4.2) 


The  radiation  impedance  for  the  monopole  sphere  which  has  a  uniform,  radial  velocity  distri- 
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bution  and  uniform  pressure  on  the  surface  is  given  as  [4] 


Zr  =  Ana2  pc 


(ka)2+ika 
1 +(ka)2 


(4.3) 


and  the  velocity  over  the  spherical  transducer  surface  can  be  calculated  by  using  an  electrical 
equivalent  circuit  for  the  spherical  transducer  (Figure  2.2).  In  this  figure  of  the  equivalent 
circuit,  the  right  side  of  the  transformer  contain  the  components  of  the  mechanical  impedance. 
Rm  is  negligible  compared  to  the  stiffness  and  the  compliance,  so  the  mechanical  impedance 
becomes  approximately 


Zm  =  i\  coM 


(DCZ 


(4.4) 


with  velocity 

NV 

u  =  - - -  , 

(Z/-  +  Zw) 


(4.5) 


where  NV  is  the  piezoelectric  driving  force. 


The  lumped  mass,  M,  the  mechanical  compliance,  C^,  and  the  electromechanical  turns  ratio,  N, 
can  be  obtained  from  equation  (2.23). 


The  total  surface  acoustic  radiated  pressure  is  calculated  by  substitution  of  equation  (4.3)  and 
(4.5)  into  equation  (4.2).  Comparison  of  the  results  for  an  analytical  calculation  and  the  COMSOL 
computation  are  shown  in  Table  4.4. 


Method 

Real 

Imag 

Amp 

approximate  analytical  solution 

-9.6 

17.88 

20.29 

COMSOL  Result 

-9.4 

17.11 

19.52 

Table  4.4:  Comparison  of  the  free-field  radiated  pressure  at  the  surface  of  the  spherical  trans¬ 
ducer  computed  in  COMSOL  against  those  calculated  from  an  equivalent  circuit. 


The  results  in  Table  4.4  show  that  the  COMSOL  result  and  the  approximate  analytical  solution 
are  in  agreement.  The  magnitude  value  of  the  COMSOL  result  differs  by  3.7  percent  from  the 
approximate  analytical  solution.  This  is  an  acceptable  difference,  since  the  approximate  ana¬ 
lytical  solution  is  derived  from  the  equivalent  circuit  model  of  a  spherical  transducer,  based  on 
the  assumption  that  its  small  cross-sectional  dimensions  are  negligible  when  compared  with  its 
radius  and  with  no  internal  losses  [20],  whereas  the  COMSOL  model  is  fully  three-dimensional. 
The  approximate  solution  does  provide  insight  as  to  the  physical  behavior  of  the  transducer. 
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4.2.6  T-Matrix  Calculation 

Nine  different  spherical  harmonic  components  (n=0,l,2  m=-n..n)  of  an  incident  wave  pressure 
were  used  in  defining  the  background  pressure  field.  The  amplitude  of  a  given  incident  wave 
pressure  mode  was  set  to  “1”  with  the  remaining  mode  amplitudes  set  to  zero  and  scattered 
component  amplitudes  computed  (see  Section  2.4.3  for  details).  For  example,  to  determine  the 
first  column  of  the  T-matrix,  an  incident  wave  with  P(')0  =  1  and  all  other  P'nm  =  0  is  imposed. 
Then  resulting  column  matrix  of  scattered  wave  component  calculated  using  (4.6) 


and 


C2  I 

D  _  '~'nmlnm  (A  ^ 

K-nm  —  (4.0) 

K  ( ka ) 


... 

Roo  . 

1 

ps 

1  nm 

Rnm  . 

0 

Evaluation  of  the  surface  integrals  were  performed  at  the  surface  of  the  spherical  thin-shell 
transducer  using  the  solution  data  set.  One  column  of  the  T-matrix  was  derived  from  each 
computation,  so  nine  computations  were  undertaken  in  order  to  complete  the  T-matrix.  The 
T-matrix  of  a  spherical  thin-shell  shell  transducer  is  given  in  Table  4.5.  The  accuracy  of  the  T- 


T-matrix  Elements 


pi 

M)0 

pi 

M-l 

pi 

MO 

pi 

Ml 

pi 

M— 2 

pi 

M-l 

pi 

MO 

pi 

Ml 

pi 

M2 

ps 

MX) 

6.91E-02 

2.60E-08 

1.18E-08 

7.36E-08 

7.51E-08 

4.51E-07 

1.94E-06 

6.99E-06 

2.81E-06 

ps 

M-l 

3.23E-08 

7.43E-02 

3.32E-08 

5.22E-08 

8.74E-08 

3.99E-08 

8.45E-07 

2.41E-06 

3.98E-06 

ps 

MO 

4.33E-09 

1.37E-08 

7.43E-02 

3.31E-08 

1.13E-07 

8.15E-08 

3.44E-07 

5.07E-07 

2.47E-06 

ps 

Ml 

9.92E-09 

2.10E-09 

1.37E-08 

7.43E-02 

1.03E-07 

1.84E-07 

3.66E-07 

9.32E-08 

1.40E-06 

M  2 

7.00E-07 

1.19E-06 

2.18E-06 

3.50E-06 

9.21E-01 

1.01E-04 

1.07E-04 

5.55E-05 

9.18E-04 

Mi 

2.61E-06 

9.96E-08 

3.85E-07 

2.15E-06 

2.58E-05 

9.21E-01 

5.22E-05 

1.14E-04 

5.55E-05 

ps 

MO 

8.93E-07 

3.21E-07 

3.05E-07 

7.39E-07 

4.20E-06 

8.07E-06 

9.21E-01 

5.22E-05 

1.07E-04 

ps 

Ml 

1.63E-07 

1.64E-07 

6.13E-08 

4.53E-08 

5.20E-07 

2.14E-06 

8.07E-06 

9.21E-01 

1.01E-04 

ps 

M2 

3.73E-08 

9.21E-08 

9.97E-08 

7.33E-08 

1.80E-06 

5.20E-07 

4.20E-06 

2.58E-05 

9.21E-01 

Table  4.5:  The  magnitude  of  T-matrix  elements  of  a  spherical  thin-shell  transducer  computed  in 
COMSOL 


matrix  can  be  determined  by  examining  the  non-diagonal  elements  because  the  applied  single¬ 
component  of  incident  wave  should  theoretically  result  in  only  the  same  harmonic  component 
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of  the  scattered  wave.  Therefore,  the  diagonal  values  of  T-matrix  should  be  non-zero  for  a 
spherical  thin-shell  transducer,  while  the  off-diagonal  values  should  vanish  [12, 14]. 

To  investigate  the  off-diagonal  elements  accurately,  the  magnitude  of  the  off-diagonal  elements 
in  the  T-matrix  were  normalized  by  the  magnitude  of  diagonal  elements  in  the  same  column 
and  the  result  is  shown  in  Table  4.6.  The  magnitude  of  the  diagonal  elements  are  significantly 
greater  than  non-diagonal  elements  as  expected.  The  largest  “leakage”  is  1.10E-04  and  occurs 


Normalized  T-matrix 


pi 

MX) 

pi 

M-l 

pi 

MO 

pi 

Ml 

pi 

M2— 2 

pi 

r2-l 

pi 

M20 

pi 

r2l 

pi 

r22 

Poo 

1 

3.50E-07 

1.58E-07 

9.92E-07 

8.15E-08 

4.90E-07 

2.10E-06 

7.59E-06 

3.05E-06 

PI- 1 

4.67E-07 

1 

4.47E-07 

7.03E-07 

9.48E-08 

4.34E-08 

9.18E-07 

2.62E-06 

4.32E-06 

ps 

MO 

6.26E-08 

1.85E-07 

1 

4.46E-07 

1.23E-07 

8.85E-08 

3.74E-07 

5.51E-07 

2.68E-06 

ps 

Ml 

1.43E-07 

2.82E-08 

1.84E-07 

1 

1.12E-07 

2.00E-07 

3.97E-07 

1.01E-07 

1.52E-06 

pi-2 

1.01E-05 

1.60E-05 

2.93E-05 

4.72E-05 

1 

1.10E-04 

1.17E-04 

6.02E-05 

9.96E-04 

ps 

M2— 1 

3.78E-05 

1.34E-06 

5.18E-06 

2.90E-05 

2.80E-05 

1 

5.66E-05 

1.24E-04 

6.02E-05 

ps 

r20 

1.29E-05 

4.33E-06 

4.11E-06 

9.95E-06 

4.56E-06 

8.76E-06 

1 

5.66E-05 

1.17E-04 

ps 

r21 

2.36E-06 

2.21E-06 

8.25E-07 

6.10E-07 

5.64E-07 

2.33E-06 

8.76E-06 

1 

1.10E-04 

ps 

r22 

5.40E-07 

1.24E-06 

1.34E-06 

9.87E-07 

1.96E-06 

5.64E-07 

4.56E-06 

2.80E-05 

1 

Table  4.6:  The  magnitude  of  the  non-diagonal  elements  in  the  T-matrix  were  normalized  to  the 
magnitude  of  diagonal  elements  in  the  same  column 


with  quadrupole-quadrupole  spherical  harmonic  component  interactions.  It  is  remarkable  that 
the  leakages  are  so  small  and  demonstrate  the  quality  of  the  T-matrix  calculations. 

The  diagonal  elements  of  the  T-matrix  are  shown  in  Table  4.7.  For  the  spherical  thin-shells. 


ELEMENT 

REAL  PART 

IMAG  PART 

MAGNITUDE 

PHASE(DEG) 

Tn 

-4.778E-03 

-6.895E-02 

6.911E-02 

-9.40E+01 

T22 

-5.514E-03 

7.405E-02 

7.426E-02 

9.43E+01 

T33 

-5.514E-03 

7.405E-02 

7.426E-02 

9.43E+01 

T44 

-5.514E-03 

7.405E-02 

7.426E-02 

9.43E+01 

T55 

-8.462E-01 

-3.642E-01 

9.212E-01 

-1.57E+02 

T<36 

-8.462E-01 

-3.641E-01 

9.212E-01 

-1.57E+02 

T11 

-8.462E-01 

-3.641E-01 

9.212E-01 

-1.57E+02 

Tss 

-8.462E-01 

-3.641E-01 

9.212E-01 

-1.57E+02 

T99 

-8.462E-01 

-3.642E-01 

9.212E-01 

-1.57E+02 

Table  4.7:  Diagonal  T-matrix  elements  for  a  spherical  thin-shell  transducer 


the  off-diagonal  elements  of  the  T-matrix  computed  with  COMSOL  were  fairly  close  to  zero. 


39 


Therefore,  the  off-diagonal  elements  were  neglected  in  array  calculations  and  the  T-matrix  for 
the  spherical  thin-shell  transducer  is  expressed  as 


TLSphere  - 


-0.005 -0.069i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.006 +  0.074i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.006 +  0.074i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.006  +  0.0741 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846- 

-0.3641  0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846  -0.3641 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846 -0.364i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846 -0.364i 

0 

0 

0 

0 

0 

0 

0 

0 

0 

-0.846 

4.3  The  Flextensional  Ring-Shell  Transducer 

The  second  major  objective  of  this  thesis  is  to  model  the  performance  of  a  linear  array  of  flex¬ 
tensional  ring-shell  transducer  elements.  As  no  analytical  solution  is  available,  an  accurate  3D 
model  must  be  constructed  to  calculate  the  T-matrix  elements  of  the  class  V  flextensional  ring- 
shell  transducer  in  COMSOL.  The  excellent  result  obtained  for  the  thin  shell  FEM  calculations 
lent  confidence  that  high  accuracy  results  could  be  obtained  for  the  ring  shell  transducer  as  well. 
Results  obtained  with  COMSOL  were  compared  to  prior  work  as  well  as  to  the  manufacturer’s 
data  in  order  to  investigate  the  accuracy  of  the  model  [7, 15]. 

Previous  research  is  provided  in  the  thesis  work  of  Rogerio  Pinto  [7].  It  should  be  noticed  that 
an  ultimate  goal  of  this  research  is  to  model  the  performance  of  the  “double-billboard“  array  of 
flextensional  transducers  that  was  built  and  tested  circa  1990  [1]. 


4.3.1  Geometry 

The  dimensional  information  of  the  flextensional  ring- shell  transducer  differs  slightly  from  that 
of  prior  research.  The  dimensional  nodes  from  the  previous  attempt  to  model  the  3D  flex¬ 
tensional  ring-shell  in  ATILA  FEM  were  used  as  a  reference  for  the  present  geometry.  The 
coordinate  values  of  the  nodes  were  imported  to  Microsoft  Excel  and  the  required  dimensions 
of  the  geometry  were  then  computed. 

The  elements  were  built  using  3D  geometric  primitives  in  COMSOL  (e.g.,  sphere,  circle,  etc.)  in 
order  to  simplify  the  geometry.  The  ring-shell  transducer  was  located  inside  of  a  fluid  spherical 
domain  with  radiating  boundaries  at  a  total  radius  of  2.5  m. 

The  dimensions  of  the  geometry  are  shown  in  Table  4.8. 
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Parameter 

Value  [cm] 

Shell  outer  radius 

67.32 

Shell  height  from  symmetry  center 

14.3 

Shell  thickness 

1.04 

Piezoelectric  ring  edge 

8.0 

Inner  radius  of  piezoelectric  ring 

31.9 

Outter  radius  of  piezoelectric  ring 

39.9 

Fiber  wrap  height 

8.0 

Fiber  wrap  width 

0.5 

Radius  of  water  domain 

250 

Table  4.8:  The  ring-shell  transducer  dimensional  parameters 


4.3.2  Material 

Piezoelectric  ceramic  ring 

The  piezoelectric  ceramic  ring  has  144  tangentially-poled  lead  zirconate  titanate  ceramic  plates, 
separated  by  72  steel  wedges.  The  modeling  for  each  tangentially-poled  ceramic  plate  requires 
coordinate  system  orientation  (Figure  3.2)  and  this  process  increases  the  computation  time. 
Therefore,  “smeared”  material  properties  were  used  for  the  composite  piezoelectric  ceramic 
ring  and  steel  wedges  in  simplifying  the  model.  Also,  the  polarization  direction  altered  from 
the  tangential  direction  into  the  axial  direction. 

The  properties  of  the  smeared  piezo-ring,  shown  in  Figure  4.7  were  obtained  by  McMahon 
and  Armstrong  through  in-air  measurements  of  the  segmented  ring  during  various  stages  of 
assembly  [7].  The  effect  of  compression  due  to  fiberglass  wrapping  and  material  losses  were 
included  in  the  piezoelectric  material  properties. 


Figure  4.7:  Smeared  piezo-ring 


41 


A  further  transformation  is  required  of  the  elasto-piezoelectric  material  property  matrix  in  order 
to  model  a  circumferentially-polarized  ring  as  an  axially-polarized  ring.  The  required  steps  are 
explained  in  Appendix  A.  After  transformation,  the  elasto-piezoelectric  matrix  becomes 
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for  which  the  smeared  material  properties  are: 


Parameter  Value  Unit  Notation 


Density 

Compliance  Matrix 

7550 

kg/m3 

sn 

8.202-0.32808i 

10“12m2/N 

**12 

-2.703+0. 108 12i 

10_12m2/N 

.E 

**13 

-3.538+0.14152i 

10_12m2/N 

**33 

11.99-0.4796i 

10_12m2/N 

**44 

26-1.04i 

10-12m2/N 

S66 

21.8-8.72i 

10-12m2/N 

Coupling  Matrix 

d3\ 

-1 

10_10C/N 

d33 

2.139 

10-10C/N 

d\5 

-2.811 

10-10C/N 

Relative  Permittivity 

efi/eo 

1297-5.1861 

4/®0 

1297-5.1861 

£o 

8.85 

10“12F/M 

Table  4.9:  The  flextensional  transducer  smeared  piezoelectric  material  properties 


To  select  an  axial  polarization  direction  the  user-defined  base  vector  system  was  added  in  the 
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definition  branch  of  the  model  and  the  cylindrical  coordinate  system  was  defined  using  the  unit 
basis  vectors  (i,  j,  k)  [21] 

xi  =  er  =  cos  (f)i  +  sin  </> j 

x2  =  =  —  sin  <p  i  +  cos  (j)j  (4.9) 

x3  =  ez  =  k  . 

It  should  be  noted  that  this  modification  of  polarization  requires  adjustment  of  the  computed 
pressure  values.  The  required  correction  was  described  in  the  masters  thesis  of  Rogerio  Pinto 
as  the  following  step  [7] 

“divide  the  displacement  and  pressure  fields  by  the  ratio  of  the  circumferential  length  of  a 
‘smeared’  piezoelectric  element  (which  is  equal  to  1/144  of  the  the  ring  circumferential  length) 
to  its  height.  In  the  present  case  this  ratio  is  0.1958.” 

The  correction  factor  calculated  using  this  ratio  as  -14.16  dB,  i.e.  a  value  of  14.6  dB  was  added 
to  the  computed  pressures,  for  a  given  voltage  drive.  Note  this  correction  is  only  required  for 
computing  the  radiated  pressure  field;  it  does  not  apply  for  the  short-circuit  scattering  calcula¬ 
tions  required  to  compute  the  T-matrix  values. 

Shell 

The  shell,  shown  in  Figure  4.8  was  made  of  a  special  alloy  steel  (ST4340)  and  modeled  by  the 
user  defined  material  properties  given  in  Table  4.10,  and  taken  from  an  input  data  file  of  ATILA 
FEM  [7], 


Figure  4.8:  ST4340  shell 


Fiberglass 

The  fiberglass  sheath,  shown  in  Figure  4.9,  provides  compression  to  the  ring,  having  been  ap¬ 
plied  under  tension  and  bonded  with  epoxy  resin,  and  the  smeared  piezoelectric  ceramic  proper- 
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Parameter 

Value 

Notation 

Density 

7700  kg/m3 

P 

Young’s  Modulus 

195  GPa 

Y 

Poisson’s  ratio 

0.28 

V 

Table  4.10:  Shell  material  properties 


ties  include  the  effect  of  the  fiberglass  wrapping  compression.  For  that  reason,  fiberglass  wrap¬ 
ping  is  added  without  any  specified  compression  to  the  model;  Only  its  linear  elasto-dynamic 
properties  are  specified,  as  shown  in  Table  4.11. 


Parameter 

Value 

Notation 

Density 

2000  kg/m3 

P 

Young’s  Modulus 

10  GPa 

Y 

Poisson’s  ratio 

0.28 

V 

Table  4.11:  Fiberglass  material  properties 
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Fluid  Domain 

The  fluid  domain  has  the  same  material  properties  as  given  in  Table  4.2.  The  fluid  media  include 
exterior  water  and  interior  air. 

4.3.3  Meshing 

The  maximum  element  size  is  defined  with  respect  to  the  Nyquist  criteria  which  states  that  the 
maximum  element  size  must  be  smaller  than  A/5  (This  criteria  was  applied  with  respect  to 
the  speed  of  sound  in  water).  The  frequency  used  is  605  Hz  which  is  close  to  the  operating 
frequency  and  also  provides  “ ka=V\ 

The  upper  surface  of  the  piezo-ring  was  meshed  with  free  triangular  mesh  elements  using  the 
size  adjustment  of  “extra  fine.”  The  meshed  surface  was  swept  by  a  built-in  feature  of  COMSOL 
through  the  axial  direction  to  the  bottom  surface  of  the  piezo-ring. 


(a)  (b) 

Figure  4.10:  (a)  Converted  faces  with  quadrilateral  mesh  (b)  Swept  mesh  with  free  triangular 
elements 

The  number  of  the  elements  swept  in  were  adjusted  to  be  10  equal  hexahedron  elements  using 
the  “distribution”  node  feature  of  COMSOL  (Figure  4.10b).  The  surface  boundary  between  ring 
faces  and  the  surrounding  fluid  was  converted  to  a  quadrilateral  mesh  to  generate  an  adjacent- 
free  tetrahedral  mesh  for  the  surrounding  domain  (Figure  4.10a). 
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The  remaining  geometry  was  meshed  using  free  triangular  mesh  elements  with  a  size  adjustment 
of  “finer”  (Figure  4.1 1).  The  complete  mesh,  shown  in  Figure  4.12  consists  of  207,003  elements 
and  the  average  element  quality  of  the  entire  geometry  is  0.66. 


Figure  4.12:  Complete  mesh  of  entire  geometry 


Geometry  Entity  Level 

Element  Type 

Amount 

Avg.  Element  Quality 

All  Domain 

Tetrahedral 

200,483 

0.67 

Pyramid 

2,680 

0.20 

Prism 

3,840 

0.48 

Boundaries 

Radiation 

Triangular 

5,056 

0.99 

Piezo  ceramic  ring 

Triangular 

3,496 

0.54 

Wet  surface  of  transducer 

Triangular 

7,588 

0.91 

Entire  Geometry 

207,003 

0.66 

Table  4.12:  Meshing  properties  of  the  3D  spherical  transducer  model  with  extra  fine  mesh  for 
the  piezo-ring  and  with  a  size  adjustment  of  finer  for  the  remaining  geometry 


4.3.4  Validation  of  Model 

To  quantify  validation  of  the  model,  the  source  level  (SL)  at  3000  Vrms  at  1  m  was  computed 
in  COMSOL,  and  the  correction  factor,  -14.16  dB,  was  added  to  the  computed  SL  (See  Chapter 
4.3.2  for  details). 
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Figure  4.13  shows  the  SL  as  a  function  of  the  frequency  computed  by  COMSOL  using  the 
baseline  properties  at  10  Hz  intervals,  from  200  Hz  to  2000  Hz  and  Figure  4.14  shows  a  com¬ 
parison  of  the  SL  computed  in  ATILA  FEM  [7]  against  the  manufacturer’s  data  of  the  SL  as  a 
function  of  frequency.  The  number  of  DOFs  solved  was  495,594  and  the  solution  time  for  each 
frequency  was  255  s. 

The  first  peak  value  of  21 1.6  dB  re  1  /iPa/V  at  1  m  occurred  at  605  Hz  with  3000  Vrms  driven 
voltage.  While  this  value  differs  by  1.4  dB  in  SL  from  the  manufacturer’s  data  (213  dB),  it  is 
much  better  than  ATILA  FEM  results  which  differed  by  4  dB  from  the  measured  value  at  the 
resonance  frequency. 


SOURCE  LEVEL  (SL)  DRIVEN  AT  3000Vrms 


Figure  4.13:  SL  at  lm  on  acoustic  axis  driven  3000  Vrms,  obtained  with  COMSOL  ACPZ 
interface  (main  model). 


The  very  narrowband  irregularities  observed  in  Figure  4.13  (660  Hz,  ~1200  Hz,  ~  1700  Hz  ) 
were  found  to  be  uncontrolled  shell  modes.  They  are  not  seen  in  experimental  measurements. 

Figure  4.15  depicts  the  bandwidth  at  resonance  frequency,  and  it  shows  the  interpolated  value 
of  ±3 dB  with  respect  to  computed  values.  The  spurious  value  at  660  Hz  was  omitted  as  being 
misleading.  The  mechanical  quality  factor,  Qm,  was  calculated  4.35  and  a  bandwidth  of  139 
Hz.  The  Qm  differ  0.4  and  the  bandwidth  16  Hz  from  the  manufacturer’s  measured  value. 

The  model  shows  (Figure  4.13)  the  second  resonance  of  the  developed  model  occurring  at  950 
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Figure  4.14:  SL  on  acoustic  axis  driven  3000  Vrms,  obtained  from  the  master  thesis  of 
Pinto,  1993  [7] 


TRANSMITTING  VOLTAGE  RESPONSE(TVR)  DRIVEN  AT  3000Vrms 


Figure  4.15:  The  bandwidth  at  primary  resonance  frequency(main  model) 


Hz  which  is  in  excellent  agreement  with  manufacturer’s  data,  but  the  model  displays  an  approx¬ 
imately  6  dB  higher  peak  after  the  secondary  resonance  frequency.  As  the  frequency  increases 
to  1000  Hz,  the  model  shows  some  discrepancies.  After  1600  Hz  the  model  shows  a  lack  of 
convergence  to  manufacturer’s  data. 
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At  first  glance,  the  expected  reason  for  this  deviation  may  be  due  to  internal  material  losses.  The 
internal  material  losses  were  not  included  in  the  model  and  could  vary  the  computed  SL  values. 
Secondly,  the  mesh  quality  (Figure  4.16)  may  affect  the  results  at  higher  frequencies  since 
meshing  quality  is  related  to  operating  wavelength.  The  proper  meshing  size  at  low  frequency 
becomes  too  coarse  at  high  frequency.  Therefore,  increasing  the  frequency,  reduces  the  mesh 
quality  of  the  model.  It  should  be  noted,  that  the  mesh  size  of  the  developed  model  was  adjusted 
according  to  “/= 474  Hz”  and  the  selected  mesh  size  was  not  refined  enough  for  high  frequency. 

Further  investigation  of  meshing  quality  indicated  another  issue:  the  meshing  of  thin  structures 
deteriorated.  Coupling  of  a  thin  structure  with  a  larger  object  obstructs  the  growth  rate  of  the 
mesh  size.  A  sample  of  this  situation  occurs  with  the  coupling  between  fiberglass  wrapping  and 
the  piezo-ring  domain.  The  quality  distribution  of  the  mesh  is  shown  in  Figure  4.16. 


0.3 

0.2 

0.1 

0 

▼  0.0155 

Figure  4.16:  Mesh  quality  plot  of  the  main  model  for  the  ring-shell  transducer. 

Although,  the  model  shows  some  deviation  above  the  secondary  resonance  frequency,  the  model 
still  demonstrates  quite  good  agreement  with  the  manufacturer’s  data  between  400-1000  Hz 
(“ka=  1”)  and  better  agreement  than  the  previous  model  (Pinto, 1993). 

The  developed  model  is  sufficient  for  the  purpose  of  a  T-matrix  calculation  because  the  perfor¬ 
mance  of  the  class  V  flextensional  transducer  will  reside  in  vicinity  of  “ka=  1”,  and  so  was  used 
as  the  base  model  for  the  T-matrix  calculation. 

4.3.5  T-Matrix  Calculation 

The  spherical  harmonic  components  of  an  incident  standing  wave  were  applied  at  the  fluid/solid 
surface  of  the  transducer,  and  the  spherical  harmonic  scattered  amplitudes  were  determined 
from  the  resultant  scattered  pressure  field  in  the  fluid. 

Post  processing  properties  were  used  to  create  a  spherical  parametrized  surface  within  the  water 
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domain  at  radii  of  0.5  m.  Figure  4.17  shows  the  defined  spherical  parametrized  surface  in  water 
domain. 


The  magnitude  of  the  T-matrix  elements  are  provided  in  Table  4.13. 


Figure  4.17:  A  parametrized  spherical  surface  was  created  at  radii  of  0.5  m  (shown  in  red;  The 
great  circles  are  on  the  radiation  boundary  at  2.5  m). 
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Table  4.13:  The  magnitude  of  T-matrix  elements  of  a  class  V  flextensional  ring-shell  transducer 
computed  in  COMSOL 


For  convenience  of  analysis,  the  magnitude  of  the  non-diagonal  elements  in  the  T-matrix  were 


50 


normalized  to  the  magnitude  of  diagonal  elements  in  the  same  column  and  the  results  are  shown 
in  Table  4.14. 


Normalized  T-matrix 
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Table  4.14:  The  magnitude  of  non-diagonal  elements  in  the  T-matrix  normalized  to  the  magni¬ 
tude  of  diagonal  element  in  the  same  column 


In  the  vicinity  of  the  primary  resonance  frequency,  605  Hz,  the  flextensional  transducer  radiates 
both  quadrupole  and  monopole  pressure  fields.  It  is  found  that  the  linear  quadrupole  radiated 
pressure  coefficient  is  about  1/10  of  the  monopole  coefficient  (  see  section  5.2  for  details)  When 
a  monopole  field  is  incident,  the  ring-shell  generates  a  quadrupole  field  and  vice  versa,  charac¬ 
teristic  of  the  resonant  vibration.  This  general  dynamic  response  of  the  transducer  can  be  seen 
from  the  T-matrix  table.  It  is  remarkable  that  the  magnitude  of  T77  (quadrupole-quadrupole)  is 
smaller  than  the  magnitude  of  T\ 7  (quadrupole-monopole)  which  are  highlighted  in  Table  4.14. 
Thus,  the  calculated  T-matrix  captures  the  interaction  between  monopole  and  quadrupole  exci¬ 
tation. 

The  T-matrix  for  the  ring-shell  transducer  is  not  diagonal  while  the  T-matrix  of  the  spherical 
thin-shell  transducer  essentially  is.  This  is  expected  behavior  as  the  geometry  of  the  flexten¬ 
sional  transducer  is  not  spherical.  It  should  be  noted,  that  the  T-matrix  is  a  kind  of  transfer 
function.  The  diagonal  elements  are  represented  in  Table  4.15. 

Similar  to  the  spherical  thin-shell  transducer,  the  T-matrix  of  the  flextensional  ring-shell  trans¬ 
ducer  is  nearly  a  diagonal  matrix  except  for  elements  of  monopole-linear  quadrupole  interac¬ 
tions  (TuJji).  Unlike  the  spherical  thin-shell  transducers,  the  monopole-linear  quadrupole 
term  must  be  included  in  the  T-matrix.  Neglecting  small  values  in  the  matrix,  the  T-matrix  of 
the  flextensional  ring-shell  becomes, 
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ELEMENT 

REAL  PART 

IMAG  PART 

MAGNITUDE 

PHASE(DEG) 

Tn 

-9.237E-01 

-1.088E-01 

9.301E-01 

-1.73E+02 

T22 

-1.538E-02 

-4.820E-02 

5.060E-02 

-1.08E+02 

T33 

-1.652E-02 

-5.435E-02 

5.680E-02 

-1.07E+02 

T44 

-2.467E-02 

-7.134E-02 

7.548E-02 

-1.09E+02 

T55 

-1.674E-03 

-4.607E-03 

4.902E-03 

-1.10E+02 

766 

-1.660E-03 

-4.129E-03 

4.450E-03 

-1.12E+02 

T11 

-3.586E-03 

-2.242E-03 

4.229E-03 

-1.48E+02 

788 

-1.660E-03 

-4.129E-03 

4.450E-03 

-1.12E+02 

Tgg 

-1.674E-03 

-4.607E-03 

4.902E-03 

-1.10E+02 

Table  4.15:  Diagonal  T-matrix  elements  for  a  class  V  flextensional  ring-shell  transducer 


TRing  shell  — 


-0.924  -0.109i  0  0  0  0  0  -0.019 -0.002i  0  0 

0  -0.015  -0.048i  0  0  0  0  0  0  0 

0  0  -0.017 -0.054i  0  0  0  0  0  0 

0  0  0  -0.025 -0.071i  0  0  0  0  0 

0  0  0  0  -0.002 -0.005i  0  0  0  0 

0  0  0  0  0  -0.002  -  0.004i  0  0  0 

-0.095  -  0.009i  0  0  0  0  0  -0.004  -  0.002i  0  0 

0  0  0  0  0  0  0  -0.002 -0.004i  0 

00000000  -0.002  -  0.005i 


4.4  Improvement  on  the  Ring-Shell  Transducer  Model  at  Higher 
Operating  Frequencies 

Although  the  accuracy  of  the  developed  model  was  verified  up  to  1000  Hz,  corresponding  to 
the  frequency  of  interest  in  this  thesis  to  calculate  the  T-matrix  elements,  the  SL  of  the  devel¬ 
oped  model  demonstrated  some  deviation  from  experimental  values  at  higher  frequencies  (Fig¬ 
ure  4.13).  To  achieve  more  accurate  results  at  higher  frequencies,  improvements  were  made 
to  the  first  developed  model.  The  first  model  developed  is  referred  to  as  the  “main  model” 
throughout  this  section. 

4.4.1  Refined  Mesh  Model 

For  frequencies  up  to  2  kHz,  the  wavelength  of  sound  in  water  is  0.745  m  (c=1490  m/s).  The 
maximum  meshing  element  size  should  be  specified  with  respect  to  A/5  at  this  frequency. 

In  this  limit,  the  maximum  element  size  of  meshed  geometry  equals  0.149  m  which  is  much 
smaller  than  the  maximum  mesh  element  size  of  the  main  model  (0.275  m).  To  obtain  an  ac¬ 
curate  result  at  higher  frequency,  the  mesh  of  the  geometry  was  refined  and  adjusted  satisfy  the 
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A/5  criteria  at  2000  Hz.  The  refined  mesh  consist  of  413,778  elements  and  733,161  DOFs. 

The  SL  of  the  refined  mesh  model  is  shown  in  Figure  4.18  and  the  SL  was  calculated  at  10  Hz 
intervals,  from  200  Hz  to  2000  Hz.  Primary  and  secondary  resonances  were  captured  at  the 
same  frequency  as  those  in  the  main  model,  as  are  the  amplitudes.  Mesh  refinement  did  not 
produce  any  significant  improvement.  Thus,  the  general  behavior  of  SL  stayed  the  same  and  the 
refined  mesh  could  not  remove  discrepancies  between  the  graph  of  computed  SL  values  and  the 
graph  of  measured  SL  values  (Figure  4.14).  It  should  also  be  noted  that  the  computation  time 
increased  due  to  the  increased  number  of  DOFs  of  the  refined  mesh. 


Figure  4.18:  SL  at  1  m  on  acoustic  axis  is  driven  at  3000  Vrms  of  the  refined  mesh  model,  the 
mesh  size  is  adjusted  with  respect  to  2000  Hz. 

4.4.2  Improved  Model 

Upon  further  investigation,  it  was  revealed  that  the  mesh  quality  was  reduced  in  the  neigh¬ 
borhood  of  the  fiberglass  wrapping  and  the  surface  of  the  shell.  The  fiberglass  wrapping  was 
re-meshed  with  triangular  elements  and  swept  through  the  axially  direction  from  top  to  bottom 
with  a  limited  maximum  mesh  element  size  to  A/30  or  0.025  m  at  2500  Hz. 

In  addition  to  mesh  improvement,  the  spherical  wave  radiation  boundary  condition  was  replaced 
with  a  spherical  perfectly  matched  layer  (PML)  to  absorb  the  incident  waves  without  reflection 
as  they  strike  the  radiation  boundary.  The  PML  surrounded  the  fluid  domain  and  had  a  0.5  m 
thickness,  reducing  the  radius  of  fluid  domain  to  2  m.  The  updated  geometry  is  represented  in 
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Figure  4.19.  The  PML  mesh  was  swept  radially  with  a  distribution  of  ten  hexahedron  elements 
through  its  thickness,  and  the  water  domain  between  PML  and  ring-shell  was  meshed  with  a 
boundary  layer  mesh.  The  boundary  layer  mesh  is  recommended  by  COMSOL  in  order  to 
enhance  the  precision  of  the  far  field  calculation  near  the  PML  [17].  The  complete  mesh  consist 
of  238,503  elements  and  792,888  DOFs. 


PML 


NG-SHELL 


Figure  4.19:  The  geometry  of  the  improved  model,  PML  surrounded  the  fluid  domain. 


Figure  4.20  shows  the  SL  of  the  improved  model  where  again  the  SL  was  calculated  at  10 
Hz  intervals,  from  200  Hz  to  2000  Hz.  The  behavior  of  the  model  changed.  The  secondary 
resonance  frequency  shifted  left  by  35  Hz  from  the  manufacturer’s  data  to  985  Hz,  and  the  SL 
became  higher  than  the  manufacturer’s  data  by  6  dB. 

The  effects  of  the  improvements  were  seen  at  frequencies  higher  than  the  secondary  resonance. 
Beyond  that  frequency,  the  SL  agrees  with  the  manufacturer’s  data  (Figure  4.14)  and  the  dif¬ 
ference  was  reduced  to  2  dB,  showing  fairly  good  agreement  with  the  manufacturer’s  data  up 
to  1600  Hz.  For  frequencies  higher  than  1600  Hz,  the  SL  again  deviates  slightly  from  the 
manufacturer’s  data. 
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SOURCE  LEVEL  (SL)  DRIVEN  AT  3000Vrms 


Figure  4.20:  SL  on  acoustic  axis  driven  at  3000  Vrms  of  the  improved  model,  PML  and  bound¬ 
ary  layer  mesh  added. 


4.4.3  Result 

In  addition  to  the  main  model,  two  different  models  were  examined  for  the  purpose  of  improving 
accuracy  at  higher  frequencies.  Figure  4.21  depicts  the  comparison  of  SL  of  the  main  model, 
the  refined  mesh  model  and  the  improved  model. 

The  main  model  and  refined  mesh  models  performed  identically  from  400  Hz  to  800  Hz.  The 
improved  model  produced  a  lower  SL  than  the  other  models  at  higher  frequencies,  and  it  con¬ 
verged  with  the  manufacturer’s  data.  Thus,  the  improved  model  performed  better  when  com¬ 
pared  to  the  other  model,  although  it  captured  the  secondary  resonance  with  a  35  Hz  difference 
from  experiment. 

Reasons  for  the  improvement  in  the  SL  value  are  due  in  part  to  the  changes  in  the  mesh  structure 
and  applied  PML  as  boundary  condition. 


SOURCE  LEVEL  (SL)  DRIVEN  AT  3000Vrms 


Figure  4.21:  The  SL  comparison  of  three  developed  model 
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The  detailed  SL  of  frequency  interest  in  the  neighborhood  of  the  first  resonance  is  illustrated 
in  Figures  4.23  through  4.24.  The  comparison  of  the  bandwidth  and  the  Q  values  are  shown  in 
Table  4.16.  The  three  models  performed  nearly  the  same  and  show  fairly  good  agreement  with 
the  manufacturer’s  data  between  250-1000  Hz. 


Model 

BW  (Hz) 

BW  Error(Hz) 

BW(%) 

Q 

Q  Error 

Q  Error(%) 

Main 

138.5 

16.5 

10.64 

4.33 

0.43 

11 

Refined 

139.5 

15.5 

10.00 

4.30 

0.4 

10.26 

Improved 

139 

16 

10.03 

4.35 

0.45 

10.35 

Table  4.16:  comparison  of  bandwidth  and  mechanical  Q  value  from  three  model  against  manu¬ 
facturer’s  data 


All  the  models  reproduced  equally  well  the  SL  up  to  1000  Hz.  For  this  reason,  the  main  model 
was  preferred  for  the  T-matrix  element  calculation  because  of  reduced  computation  time  for 
equivalent  accuracy. 


Model 

number  of  element 

average 

worst 

DOFs 

Computation  Time 

main 

207,003 

0.66 

2.29E-05 

495594 

258  s. 

refined 

413,778 

0.73 

4.23E-05 

733161 

353  s. 

improved 

238,503 

0.56 

2.23E-06 

792888 

550  s. 

Table  4.17:  The  computation  time  values  and  the  general  mesh  information  of  the  three  model 


SL  DRIVEN  AT  3000Vrms(MAIN  MODEL) 


Figure  4.22:  The  bandwidth  at  primary  resonance  frequency  (main  model) 
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Figure  4.23:  The  bandwidth  at  primary  resonance  frequency  (refined  model) 


SL  DRIVEN  AT  3000Vrms  (IMPROVED  MODEL) 


Figure  4.24:  The  bandwidth  at  primary  resonance  frequency  (improved  model) 
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CHAPTER  5: 

Linear  Array  Performance  Modeling  by  Using  T-Matrix 

Method 


The  vertical  linear  array  with  the  spherical  thin-shell  transducer  elements  and  the  vertical  linear 
array  with  the  ring-shell  transducer  elements  were  programmed  in  MATLAB,  using  the  T-matrix 
method.  The  required  T-matrix  elements  for  both  types  of  the  transducers  were  reported  in 
Chapter  4. 

5.1  The  Development  Sequence  of  an  Array  Model 

1.  Computation  of  the  free-held  radiated  pressure  coefficients 

2.  Calculation  of  the  translation  matrix  using  the  addition  theorem 

3.  Implementation  of  the  T-matrix  to  calculate  the  scattered  pressure  coefficients 

4.  Computation  of  the  total  acoustic  pressure  held  and  generation  of  the  far-held  beam  pat¬ 
tern 


5.2  Computation  of  the  Free-Field  Radiated  Pressure  Coefficients 


The  free-held  radiated  pressure  coefficients  of  both  transducer  types  were  calculated  at  a  1  V 
peak  driven  potential.  That  a  driving  peak  was  used  versus  1  Vrms  is  of  no  consequence  as 
only  normalized  beam  patterns  are  presented  in  this  chapter.  The  free-held  radiated  pressure 
amplitudes  were  computed  by  solving  equation  (2.40)  using  COMSOL.  The  number  of  spher¬ 
ical  harmonic  components  kept  was  truncated  at  the  the  quadrupolar  degree,  n- 2,  so  the  total 
number  of  spherical  harmonic  components  kept  was  9  (n=0, 1 ,2  and  m=-n..n). 

Results  for  the  spherical  thin-shell  transducer  are  listed  in  Table  5.1,  and  the  results  for  the 
ring-shell  transducer  are  given  in  Table  5.2. 
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operation  freq 

474Hz(ka=l) 

Mode 

Real 

Imaginary 

Amplitude 

Amplitude 

lpool 

Angle(rad) 

PV 

MO 

1.349E+00 

1.946E+01 

1.951E+01 

1.000E+00 

1.502E+00 

pr 

M-l 

6.323E-06 

6.863E-06 

9.331E-06 

4.784E-07 

8.263E-01 

pr 

MO 

3.608E-07 

1.330E-06 

1.378E-06 

7.067E-08 

1.306E+00 

pr 

Ml 

2.909E-06 

-9.429E-07 

3.058E-06 

1.568E-07 

-3.135E-01 

pr 

M-2 

2.633E-04 

1.202E-04 

2.894E-04 

1.484E-05 

4.285E-01 

pr 

M-l 

-6.943E-04 

3.410E-04 

7.735E-04 

3.966E-05 

2.685E+00 

pr 

MO 

2.027E-04 

9.465E-05 

2.237E-04 

1.147E-05 

4.369E-01 

pr 

Ml 

3.074E-05 

3.885E-05 

4.954E-05 

2.540E-06 

9.015E-01 

pr 

M2 

-1.306E-06 

-8.569E-06 

8.668E-06 

4.444E-07 

-1.722E+00 

Table  5.1:  The  radiated  pressure  coefficients  of  the  spherical  thin-shell  transducer  for  monopole, 
dipole,  and  quadrupole  (n- 0,  1,  2  and  m=-n..n),  for  1  V  peak  drive. 


operation  freq 

60 5 Hz  (ka=l) 

Mode 

Real 

Imaginary 

Amplitude 

Amplitude 

lpool 

Angle  (rad) 

pr 

MO 

6.431E+00 

9.823E-01 

6.506E+00 

1.000E+00 

1.516E-01 

pr 

M-l 

1.026E-04 

-6.191E-04 

6.275E-04 

9.645E-05 

-1.407E+00 

pr 

MO 

1.109E-03 

3.419E-03 

3.594E-03 

5.524E-04 

1.257E+00 

pr 

Ml 

4.475E-04 

-4.552E-04 

6.384E-04 

9.812E-05 

-7.939E-01 

pr 

M— 2 

-4.804E-06 

1.083E-03 

1.083E-03 

1.665E-04 

1.575E+00 

pr 

M-l 

1.188E-04 

1.791E-04 

2.149E-04 

3.303E-05 

9.853E-01 

pr 

MO 

6.738E-01 

3.858E-02 

6.749E-01 

1.037E-01 

5.719E-02 

pr 

Ml 

-2.213E-05 

2.659E-05 

3.459E-05 

5.317E-06 

2.265E+00 

pr 

M2 

8.583E-06 

-6.395E-06 

1.070E-05 

1.645E-06 

-6.403E-01 

Table  5.2:  The  radiated  pressure  coefficients  of  the  ring-shell  transducer  for  monopole,  dipole, 
and  quadrupole  (n- 0,  1,  2  and  m=-n..ri),  for  1  V  peak  drive. 


Table  5.1  indicates  that  the  free-field  radiated  pressure  amplitudes  of  the  spherical  thin-shell 
transducer  are  negligible  except  for  the  monopole  (7q0),  as  expected.  Because  of  symmetry, 
coefficients  other  than  P[)0  are  close  to  zero,  indicating  the  accuracy  of  the  finite  element  model. 

Similarly,  Table  5.2  shows  that  the  free-field  radiated  pressure  amplitudes  of  the  ring-shell  trans¬ 
ducer  are  negligible  except  for  the  monopole  (Pq0)  and  linear  quadrupole  {P^ 0),  the  monopole 
coefficient  being  about  “10”  times  the  linear  quadrupole  coefficient.  This  is  expected  due  to  the 
nature  of  the  flexural  motion  of  the  ring- shell  transducer. 
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Negligible  values  of  the  free-field  radiation  are  ignored  in  far-field  calculation.  The  free-field 
radiated  pressure  vectors  then  become  for  a  1  V  peak  drive, 


pR  — 

1  Sphere 


1.345 +19.459i 
0 
0 
0 
0 
0 
0 
0 
0 


PR 

Ring— shell 


6.431+0.982 i 
0 
0 
0 
0 
0 

0.674  +  0.039/ 
0 
0 


and  they  were  imported  to  MATLAB  for  further  calculation. 


(5.1) 


(5.2) 


5.3  Calculation  of  the  Translation  Matrix 

The  translation  matrix, “Giq”,  is  needed  to  relate  the  outgoing  wave  component  from 

transducer  “q”  to  the  “n,  mth”  free-field  incident  component  on  transducer  The  translation 
matrix  is  obtained  from  a  spherical  addition  theorem  (See  Chapter  2.4.5  for  details).  Bold 
capital  “G”  denotes  a  translation  matrix  from  one  transducer  to  another  in  the  array 
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Gz,  oo  Gi,  oo  Gz,oo  Gi,  oo  Gz,oo  Gz,oo  G/,oo  G/,oo  Gz,oo 

4,00  4,1  —  1  9,10  9,11  9,2—2  9,2—1  4,20  9,21  9,22 

G/,i-i  G/,i-i  G/,i-i  G/,i-i  G/,i-i  G/,i-i  G/,i-i  Gz,i-i  G/,1-1 

9,00  9,1-1  <7,10  9,11  q,  2—2  4, 2-1  4, 20  9, 21  9,22 

G/,10  G/,10  G/,10  G/,10  Gz, 10  Gz, 10  Gz, 10  Gz, 10  Gz, 10 

9,00  9,1-1  4,10  9,11  q,  2-2  4,2-1  4,20  4, 21  4,22 

Gz,n  Gz,n  Gz,n  Gz,n  Gz,n  Gz,n  Gz,n  Gz,n  Gz,n 

4,00  9,1  —  1  4,10  4,11  9,2—2  9,2—1  9,20  9,21  9, 22 

Gz,  2-2  Gz,  2-2  Gz,  2-2  Gz,  2-2  Gz,  2-2  Gz,  2-2  Gz,  2-2  Gz,  2-2  Gz,  2-2 

9,00  9,1  —  1  4,10  4,11  9,2—2  9,2—1  4,20  4,21  4,22 

Gz, 2-1  Gz, 2-1  Gz, 2-1  Gz, 2-1  Gz, 2-1  Gz, 2-1  Gz, 2-1  Gz, 2-1  Gz, 2-1 

4,00  9,1  —  1  4,10  4,11  9,2—2  9,2—1  4,20  4,21  4,22 

Gz, 20  Gz, 20  Gz, 20  Gz, 20  Gz, 20  Gz, 20  Gz, 20  Gz, 20  Gz, 20 

4,00  4,1-1  4,10  4,11  4,2-2  4,2-1  4,20  4,21  4,22 

Gz, 21  Gz, 21  Gz, 21  Gz, 21  Gz, 21  Gz, 21  Gz, 21  Gz, 21  Gz, 21 

4,00  4,1-1  4,10  4,11  4,2-2  4,2-1  4,20  4,21  4,22 

G  Z,22  G  Z,22  G  Z,22  Gz,22  Gz,22  Gz,22  Gz,22  Gz,22  Gz,22 

4,00  9,1  —  1  4,10  4,11  9,2—2  9,2—1  4,20  4,21  4,22 


Equation  (2.43)  was  employed  in  MATLAB  to  construct  the  translation  matrix,  Giq,  /  7^  q  . 


The  llh  row  refers  to  the  interaction  of  the  “Zf/i”  transducer  with  the  other  transducers  in  the  array. 
Computation  time  for  producing  the  translation  matrix  depends  on  the  number  of  elements 
within  an  array 

0  G12  •  •  •  Gin 

G21  0  •  •  •  G2N 


GnI  Gn2  0 


where 


pwrt  =  pR  +  pS 


5.4  Implementation  of  the  T-matrix 

The  scattered  pressure  amplitudes  can  be  obtained  solving  (2.49).  The  product  of  both  sides 
with  the  inverse  matrix  of  (I  —  a)  gives  the  scattered  pressure  vectors. 


where  a  is  defined  for  an  array  which  is  composed  of  the  identical  elements, 


(5.4) 
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a  = 


/  0  TG12  •  •  •  TGin  \ 

TG21  0  tg2N 

\  TGni  TGn2  •••  0  / 

Considering  a  three-element  transducer  array  configuration,  for  example,  the  free-field  scattered 
pressure  amplitudes  become 


P1 

\ 

I 

-tg12 

-tg13  ' 

-1 

0 

tg12 

tg13 

/  pf 

\ 

pf 

— 

-TG21 

1 

-TG23 

TG21 

0 

TG23 

pf 

V  Pf 

J 

-TG31 

-tg32 

I 

.  TG31 

TG32 

0 

\pf 

J 

In  Chapter  4.2.6  the  complex  coefficients  of  the  T-matrix  for  the  spherical  thin-shell  transducer 
were  calculated  as 


Tsphere  — 


-0.005  -0.069i  0  0  0  0 

0  -0.006 +  0.074i  0  0  0 

0  0  -0.006 +  0.074i  0  0 

0  0  0  -0.006  +  0.074i  0 

0  0  0  0  -0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 

0  0  0  0  0 


0  0  0  0 

0  0  0  0 

0  0  0  0 

0  0  0  0 

.846  —  0.364i  0  0  0  0 

-0.846 -0.364i  0  0  0 

0  -0.846 -0.364i  0  0 

0  0  -0.846 -0.364i  0 

0  0  0  -0.846 -0.364i 


In  Chapter  4.3.5  the  complex  coefficients  of  T-matrix  for  the  ring-shell  transducer  were  calcu¬ 
lated  as 


TRing  shell  — 


-0.924 -0.109i  0  0  0  0  0  -0.019  -  0.002i  0  0 

0  -0.015  -0.048i  0  0  0  0  0  0  0 

0  0  -0.017 -0.054i  0  0  0  0  0  0 

0  0  0  -0.025 -0.071i  0  0  0  0  0 

0  0  0  0  -0.002 -0.005i  0  0  0  0 

0  0  0  0  0  -0.002  -  0.004i  0  0  0 

-0.095  -  0.009i  0  0  0  0  0  -0.004  -  0.002i  0  0 

0  0  0  0  0  0  0  -0.002 -0.004i  0 

00000000  -0.002  -  0.005i 
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5.5  Computation  of  the  Total  Acoustic  Pressure  Field  and 
Far-Field  Beam  Pattern 

In  the  array  models,  a  line  of  N\  transducer  sources  with  the  adjacent  elements  equally  spaced 
were  analyzed  and  are  as  shown  in  Figure  5.1.  The  transducers  are  aligned  with  the  polar  axis 
(z  axis).  The  beam  patterns  of  the  developed  models,  therefore,  have  rotational  symmetry  with 
respect  to  the  array  (polar)  axis.  Beam  patterns  are  represented  by  polar  coordinates  (0,0), 
and  the  far-held  approximation  was  used  with  phase  difference  between  elements,  due  to  their 
position.  This  was  taken  into  account  for  the  calculations  of  far-held  acoustic  held  and  beam 
pattern. 


(a)  (b) 

Figure  5.1:  Orientation  of  the  transducers;  spacing  X/2  and  "ka  =  1"  (a)  Spherical  thin-shell 
transducer  (b)  Class  V  ring-shell  transducer 


The  total  acoustic  pressure  in  the  far-held  was  obtained  from  superposition  of  the  radiated  and 
the  scattered  complex  acoustic  pressure  amplitudes  originating  from  the  transducers  in  the  array 
(See  Chapter  2.4.7  for  details) .  Equations  (5.5)  and  (5.7)  were  employed  in  MATLAB  to  obtain 
far-held  pressure  amplitudes.The  radiated  pressure  amplitudes  are  reported  in  Section  5.2,  and 
the  scattered  pressure  amplitudes  can  be  obtained  from  equation  (5.4).  All  equations  were 
combined  in  MATLAB,  and  the  code  used  is  presented  in  Appendix  B. 

The  far-held  pressure  is 
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„—ikr  Nx  _  Wmax  n 

e  r'1  ,»-r0/  Y1  •(«+!) 


Pff  (r,0,4>)  =  —£<?*<»  Y,*<n+1)  Z  PL  +  PL  W(W),  (5-5) 


;=1  n=0  m=—n 


where 


k-foi—k ( xqi  sin  6 cos </>  +yoi  sin  6  sin <j>  +  zoi cos  6)  ■ 


(5.6) 


and  xqi,  you  zoi  are  the  rectangular  coordinates  of  the  origin  of  transducer 
The  normalized  far-field  beam  pattern  is 


b  =  201og10 


f  faMj )j 

\v|p//(r,0,0)|maxyl  ' 


(5.7) 


To  understand  the  array  interaction  effect  on  the  far-field  beam  pattern  different  linear  array 
configurations  were  examined.  For  all  configurations,  3D  normalized  far-field  beam  patterns 
were  plotted  using  the  T-matrix  method.  As  the  horizontal  (i.e.,  perpendicular  array  axis)  beam 
patterns  of  all  configurations  were  omni-directional,  the  detailed  investigation  of  the  effect  of 
the  T-matrix  method  on  beam  pattern  was  focused  on  the  vertical  beam  (plane  including  array 
axis  direction)  patterns.  The  resultant  normalized  vertical  beam  patterns  were  compared  against 
those  computed  for  non-interacting  simple  sources. 


5.6  Figures  of  the  Far  Field  Beam  Patterns 

The  2D  normalized  vertical  far-field  beam  patterns  were  plotted  in  dB  as  a  function  of  the  an¬ 
gle  (in  degrees).  The  beam  patterns  obtained  with  the  T-matrix  method  are  compared  to  point 
sources  on  the  same  plot.  For  the  polar  plot,  magnitudes  lower  than  “-30  dB”  were  omit¬ 
ted  since  values  lower  than  -30  dB  are  negligible  compared  to  the  maximum  intensity.  The 
“0  dB”  reference  value  was  shifted  to  30  dB  to  improve  visualization.  Similarly,  for  improving 
the  visualization  of  cartesian  plots  of  the  beam  patterns,  magnitudes  which  were  lower  than 
“-50  dB”  were  omitted  and  the  “0  dB”  reference  shifted  to  50  dB.  The  angles  “90”  and  “270” 
correspond  to  the  broadside  direction  while  “0”  and  “180”  correspond  to  the  “array  axis  (end 
fire)”  direction. 

Several  configurations  were  examined,  to  investigate  the  performance  of  the  array  and  to  deter¬ 
mine  array  interaction. 
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5.6.1  The  Model  of  the  Linear  Array  With  Spherical  Thin-Shell  Elements 

Three-Element  Array  Configuration 

Figure  5.2  depicts  the  3D  normalized  far-field  beam  pattern  of  a  linear  array  with  three  evenly- 
spaced,  identical,  spherical  thin-shell  transducer  elements  at  “ka=  1”  (/=474  Hz).  The 
inter-element  spacing  is  equal  to  A/2  where  A  is  the  acoustic  wavelength.  In  terms  of  am¬ 
plitudes,  the  broadside  is  the  major  lobe,  and  a  nodal  surface  does  not  exist.  Instead  of  the 
expected  conical  nodal  surface,  some  degradation  is  observed  between  broadside  and  end  fire. 

Figure  5.3  is  the  normalized  vertical  beam  pattern  of  the  three  evenly-spaced,  identical,  spher¬ 
ical  thin-shell  transducers.  The  inter-element  spacing  is  equal  to  A/2.  Figure  5.3b  and  Fig¬ 
ure  5.3c  compares  the  patterns  obtained  from  the  T-matrix  method  versus  the  point  sources. 
The  T-matrix  beam  pattern  shows  the  severe  interaction  at  end  fire.  The  end  fire  lobe  has  3 
dB  greater  amplitude  than  the  broadside  lobe  (main  lobe).  The  degradation  of  amplitudes  was 
observed  at  —132°,  —47°,  47°,  and  132°,  but  the  point  sources  beam  pattern  goes  to  adjacent 
nulls  at  the  same  angles  (Figure  5.3c). 

It  was  suspected  that  a  linear  quadrupolar  resonance  was  responsible  that  we  had  not  previously 
identified.  Examination  of  the  total  scattering  amplitude  revealed  very  strong  linear  quadrupolar 
scattering,  about  forty  times  as  great  as  monopolar  scattering.  An  eigenfrequency  analysis  in  the 
vicinity  of  the  operating  frequency  (474  Hz)  was  conducted  using  the  COMSOL  Multiphysics 
finite-element  code.  A  group  of  five  and  only  five  eigenfrequencies  were  found  very  close  to  474 
Hz,  at  about  477  Hz.  These  could  represent  modes  for  which  n  =  2,  m  =  -2  to  2.  Mode  shapes  for 
these  solutions  were  examined,  and  all  appeared  to  represent  quadrupole  modes,  with  four  of  the 
five  appearing  to  display  predominantly  motion  in  a  plane,  with  quadrupolar  symmetry,  and  with 
exactly  one  appearing  to  display  linear  quadrupolar  motion.  Confirmation  that  these  most  likely 
were  indeed  quadrupolar  resonances  was  obtained  by  a  second,  approximate  computation  of  the 
resonance  frequencies  of  the  axially-symmetric  eigenfrequencies  of  a  thin  elastic  shell  in  water 
[26,  pp.  281-284],  substituting  for  the  Young’s  modulus  and  Poisson’s  ratio  the  piezoelectric 
elastic  compliance  modulus  values  .vf,  and  ,vf2A>'f| .  A  linear  quadrupole  eigenfrequency  value 
of  473  Hz  was  found. 

Seven-Element  Array  Configuration 

Figure  5.4  depicts  the  3D  normalized  far-field  beam  pattern  of  the  linear  array  with  seven 
evenly-spaced,  identical,  spherical  thin-shell  transducer  elements  at  “ka=  1”  (/=474  Hz).  The 
inter-element  spacing  is  equal  to  A/2.  The  width  of  the  main  lobe  is  reduced,  and  the  directiv- 
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ity  is  enhanced  in  this  configuration  compared  to  the  three  element  array.  Two  conical  nodal 
surfaces  are  seen  between  broad  side  and  end  fire. 

Figure  5.5  is  the  normalized  vertical  beam  pattern  of  the  seven  evenly-spaced,  identical,  spher¬ 
ical  thin-shell  transducers.  The  inter-element  spacing  is  equal  to  A/2.  Figure  5.5b  and  Fig¬ 
ure  5.5c  compare  the  patterns  obtained  from  the  T-matrix  method  versus  the  point  sources.  Near 
the  broadside,  the  behavior  of  beam  patterns  is  almost  the  same  as  that  of  the  point  sources.  The 
difference  is  increased  near  and  at  the  end  fire  direction.  Separation  between  lobes  cannot  be 
detected  easily.  Three  minor  lobes  were  smeared  by  one  wide  lobe.  While  six  adjacent  nulls 
appear  in  the  beam  pattern  of  the  point  sources,  only  one  adjacent  null  occurs  in  the  T-matrix 
method. 


Eight-Element  Array  Configuration 

Figure  5.6  depicts  the  3D  normalized  far-field  beam  pattern  of  the  linear  array  with  eight  evenly- 
spaced,  identical,  spherical  thin-shell  transducer  elements  at  “£a=  1”  (/=474  Hz).  The  inter¬ 
element  spacing  is  equal  to  A/2.  Three  conical  nodal  surfaces  were  observed. 

Figure  5.7  shows  explicitly  the  difference  between  the  two  methods  for  angles  near  end  fire. 
The  T-matrix  method  does  not  produce  a  null  to  the  end  fire  direction,  while  the  point  sources 
do.  The  Figures  5.7b  and  5.7c  show  that  the  amplitude  values  of  the  T-matrix  method  away 
from  broadside  are  slightly  higher  than  the  point  sources  in  the  neighborhood  of  end  fire. 


Three-Element  Array  Configuration  With  Different  Inter-Element  Spacing  Distance 

Figure  5.3  shows  the  severe  interaction  effect  at  the  end  fire  where  the  side  lobe  has  a  larger 
amplitude  than  the  main  lobe.  To  investigate  further  the  array  interaction  of  the  three-element 
configuration,  different  inter-element  spacing  distances  for  the  three-element  array  were  exam¬ 
ined.  Figure  5.8  illustrates  the  effect  of  an  inter-element  spacing  on  the  three-element  array 
when  “ka=l.”  The  inter-element  spacing  distances  are  selected  arbitrarily  between  A/2-A/4. 
Figures  5.8a  and  5.8b  show  that  the  A/2  and  A/2.2  inter-element  spacing  distances  produce 
larger  side  lobe  than  main  lobe.  Figures  5.8c  through  5.8i  illustrate  the  shading  of  the  end  fire 
is  reduced  with  respect  to  the  decreasing  inter-element  spacing  distance,  and  it  reaches  a  min¬ 
ima  when  spacing  distance  equals  A/3.  Although  the  two  methods  show  some  difference  in 
amplitude,  they  have  the  same  beam  pattern  behavior. 
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Array  Configuration  With  Different  Numbers  of  Elements 

To  determine  the  effect  of  the  number  of  elements  in  an  array,  several  different  configurations 
were  examined.  Figures  5.9  and  5.10  represent  the  far-field  beam  patterns  for  differing  N  (num¬ 
bers  of  elements  in  the  vertical  array).  The  elements  are  spaced  equally  and  the  inter-element 
spacing  is  fixed  at  A/2.  For  convenience,  the  configurations  with  an  even  number  of  elements 
and  those  with  an  odd  number  are  analyzed  in  separate  figures. 

Figure  5.9  shows  the  end  fire  shading  is  reduced  respectively  for  an  increasing  number  of  ele¬ 
ments  in  an  array.  When  the  number  of  elements  in  an  array  is  more  than  three,  the  nulls  which 
are  predicted  near  the  end  fire  in  the  point  sources  are  shadowed.  An  analogous  shadowing 
effect  occurs  near  the  broadside  when  the  number  of  the  elements  in  an  array  is  equal  to  1 1 
(Figure  5.9e). 

Figure  5.10  shows  the  configurations  for  even-numbered  array  elements.  The  adjacent  null 
is  produced  at  end  fire  for  all  even  number  configuration  in  point  sources  due  to  the  axis  of 
symmetry  while  the  T-matrix  method  beam  patterns  display  the  array  interaction  effect,  and  the 
beam  pattern  does  not  terminate  at  the  adjacent  nulls  at  the  end  fire. 

Figures  5.11  and  5.12  demonstrate  the  3D  normalized  far-field  beam  pattern  for  the  different 
number  configurations  corresponding  to  Figures  5.9  and  5.10.  The  directionality  of  the  array 
increases  with  respect  to  N.  The  worst  directionality  appears  in  Figure  5.11a  where  the  major 
lobe  seems  to  be  at  the  end  fire.  As  shown  in  Figure  5.11c,  Figure  5. Ilf,  and  the  Figure  5.12b 
some  of  the  estimated  conical  nodal  surfaces  disappear. 

Summary  of  Results  for  Linear  Arrays  of  Spherical  Thin-Shell  Transducers 

1.  The  observation  of  the  three-element  array  configuration  shows  that  for  the  T-matrix  re¬ 
sult,  the  end-fire  lobe  is  3  dB  greater  than  the  broadside  lobe  (5.3).  The  strong  array 
interaction  was  found  to  be  due  to  strong  linear  quadrupole  scattering.  This  was  traced  to 
an  uncontrolled  (not  able  to  be  driven)  linear  quadrupole  resonance  of  the  transducer,  at 
477  Hz.  The  array  interaction  cause  a  “shading”  in  this  case  and  reduce  the  amplitude  of 
the  major  lobe. 

2.  For  the  T-matrix  results,  the  consistent  “shadowing”  effect  is  manifested  by  a  reduction 
in  destructive  interference,  and  so  true  nulls  do  not  occur. 

3.  The  “shadowing”  effect  disappears  when  the  inter-element  spacing  distance  is  A/3  for 
the  three  element  array  (Figure  5.8f). 

4.  The  configurations  of  the  linear  arrays  with  different  N  number  of  the  elements,  are  shown 
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in  Figures  5.9  and  5.10,  which  indicate  that  the  “shading”  effect  is  reduced  with  respect 
to  the  increasing  number  of  elements  while  the  “shadowing”  effect  is  increased.  It  can  be 
concluded  that  the  "shading"  effect  due  to  scattering  is  reduced  by  increasing  the  number 
of  elements  in  an  array. 

5.  It  should  be  noted  that  the  beam  pattern  of  both  methods  is  the  same  at  broadside  when 
the  number  of  elements  in  an  array  is  more  than  three  (Figures  5.9  and  5.10). 

6.  The  least  desirable  array  configuration  of  the  spherical  thin-shells  was  observed  for  a 
three-element  array  with  A/2  inter  element  spacing,  because  the  side  lobe  dominated  the 
broadside  major  lobe  (Figure  5.8a). 
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3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=3  and  d=>  /2 


1 50 


Figure  5.2:  3D  front-left  view  of  normalized  far-field  beam  pattern  of  vertical  linear  array  with 
three  identical  spherical  thin-shell  transducer  elements  obtained  with  the  T-matrix  method, inter¬ 
element  spacing  d=X/ 2,  /=474  Hz  and  ka=  1. 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=W2 


(a)  (b)  (c) 

Figure  5.3:  Normalized  far-field  vertical  beam  pattern  for  the  modeled  linear  array  with  three 
identical  spherical  thin-shell  transducer  elements  lying  along  polar  axis  (0  =  0),  inter-element 
spacing  d=X/2,  f= 474  Hz  and  ka=  1.  (a)  Orientation  of  the  spherical  thin-shell  transducers  (b) 
A  polar  plot  of  the  normalized  magnitudes  of  the  vertical  far-field  beam  pattern  (c)  A  plot  of  the 
normalized  magnitudes  of  the  vertical  far-field  beam  pattern 
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3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=7  and  d=>  12 


Figure  5.4:  3D  front-left  view  of  normalized  far-field  beam  pattern  of  vertical  linear  array  with 
seven  identical  spherical  thin- shell  transducer  elements  obtained  with  the  T-matrix  method, 
inter-element  spacing  d=X/ 2,  /=474  Hz  and  ka=  1 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=7  d=M2 


(a)  (b)  (c) 

Figure  5.5:  Normalized  far-field  vertical  beam  pattern  for  the  modeled  linear  array  with  seven 
identical  spherical  thin-shell  transducer  elements  lying  along  polar  axis  (0  =  0),  inter-element 
spacing  d=X/2,  f= 474  Hz  and  ka=  1.  (a)  Orientation  of  the  spherical  thin-shell  transducers  (b) 
A  polar  plot  of  the  normalized  magnitudes  of  the  vertical  far-field  beam  pattern  (c)  A  plot  of  the 
normalized  magnitudes  of  the  vertical  far-field  beam  pattern 
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Figure  5.6:  3D  front-left  view  of  normalized  far-field  beam  pattern  of  vertical  linear  array 
with  eight  identical  spherical  thin-shell  transducer  elements  obtained  with  the  T-matrix  method, 
inter-element  spacing  d=X/ 2,  /=474  Hz  and  ka=  1 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=8  d=>/2 


(a)  (b)  (c) 

Figure  5.7:  Normalized  far-field  vertical  beam  pattern  for  the  modeled  linear  array  with  three 
identical  spherical  thin-shell  transducer  elements  lying  along  polar  axis  (0  =  0),  inter-element 
spacing  d=X/ 2,  /=474  Hz  and  ka=  1.  (a)  Orientation  of  the  spherical  thin-shell  transducers  (b) 
A  polar  plot  of  the  normalized  magnitudes  of  the  vertical  far-field  beam  pattern  (c)  A  plot  of  the 
normalized  magnitudes  of  the  vertical  far-field  beam  pattern 
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NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/2.4 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/2.6 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/2.8 


(d)  d=X/ 2.6 


(e)  d=X/ 2.8 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=W3 


(f)  d=X/  3 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/3.2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/3.4  NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/4 


(g)  d=X  /3 .2  (h)  d=X /3 .4  (i)  d=X /4 

Figure  5.8:  The  comparison  of  the  normalized  far-held  vertical  beam  patterns  of  a  modeled 
evenly  spaced  three-element  linear  array  of  spherical  thin-shell  transducer  at  /=474  Hz  and 
ka=  1  with  different  inter-element  spacing  distance  (a-i)  Comparison  of  polar  plot  of  the  nor¬ 
malized  vertical  far-held  beam  pattern  of  T-matrix  method  against  polar  plot  of  the  normalized 
vertical  far-held  beam  pattern  of  the  point  sources. 
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NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=3  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=5  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=7  and  d=>/2 


(a)  N= 3  (b)  N=5  (c)  N=1 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=9  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=1 1  and  d=M2 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=13  and  d=>/2 


(d)  N=9 


(e)  N=\  1 


(f)  N=  13 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=15  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=17  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=19  and  d=>/2 


(g)N=l5  (h)  N=\l  (i)  A7=19 


Figure  5.9:  The  vertical  far-held  beam  patterns  for  a  modeled  vertical  linear  array  composed 
of  an  odd  number  N  of  evenly  spaced  spherical  thin-shell  transducer  elements  lying  along  the 
polar  axis  (0  =  0)  at  A/2,  /= 474  Hz  and  ka=  1  (a-i)  Comparison  of  cartesian  plot  of  the  normal¬ 
ized  vertical  far-held  beam  pattern  of  T-matrix  method  against  cartesian  plot  of  the  normalized 
vertical  far-held  beam  pattern  of  the  point  sources. 
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NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=4  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=6  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=8  and  d=>/2 


(a)  N= 4 


(b)  N= 6 


(c)  iV=8 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=10  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=12  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=14  and  d=M2 


(d)  N=  10 


(e)  N=  12 


(f)  N=  14 


(g)  ^=16 


(h)  iV=18 


(i)  Af=20 


Figure  5.10:  The  vertical  far-field  beam  patterns  for  a  modeled  vertical  linear  array  composed 
of  an  even  number  N  of  evenly  spaced  spherical  thin-shell  transducer  elements  lying  along  the 
polar  axis  (0  =  0)  at  A/2,  /'=474  Hz  and  ka=  1  (a-i)  Comparison  of  cartesian  plot  of  the  normal¬ 
ized  vertical  far-field  beam  pattern  of  T-matrix  method  against  cartesian  plot  of  the  normalized 
vertical  far-field  beam  pattern  of  the  point  sources. 
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Figure  5.11:  The  3D  normalized  far-held  beam  patterns  for  a  modeled  vertical  linear  array 
composed  of  an  odd  number  N  of  evenly  spaced  spherical  thin-shell  transducer  elements  lying 
along  the  polar  axis  (0  =  0)  at  A/2,  /= 474  Hz  and  ka= 1  (a-i)  3D  front-left  view  of  far-held 
beam  pattern  of  vertical  linear  array 
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Figure  5.12:  The  3D  normalized  far-held  beam  patterns  for  a  modeled  vertical  linear  array 
composed  of  an  even  number  N  of  evenly  spaced  spherical  thin-shell  transducer  elements  lying 
along  the  polar  axis  (0  =  0)  at  A/2,  /= 474  Hz  and  ka= 1  (a-i)  3D  front-left  view  of  far-held 
beam  pattern  of  vertical  linear  array. 
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5.6.2  The  Linear  Array  with  the  Ring-Shell  Elements 

Three-Element  Array  Configuration 

Figure  5.13  depicts  the  3D  normalized  far-field  beam  pattern  of  a  linear  array  with  three  evenly- 
spaced,  identical,  ring-shell  transducer  elements  at  uka=V ’  (/= 605  Hz).  The  inter-element  spac¬ 
ing  is  equal  to  A/2.  It  is  remarkable  that  the  conical  nodal  surfaces  do  not  degrade  relative  to 
broadside  at  end  fire 

Figure  5.14  is  the  normalized  vertical  beam  pattern  of  the  three  evenly-spaced,  identical,  ring- 
shell  transducers.  The  inter-element  spacing  is  again  equal  to  A/2.  Figures  5.14b  and  5.14c 
compare  the  beam  patterns  obtained  from  the  T-matrix  method  versus  the  point  source  method. 
Unlike  for  the  spherical  thin  shell  transducers,  the  beam  patterns  for  the  ring  shell  transducers 
are  nearly  the  same  for  both  techniques  except  true  nulls  do  not  occur  in  T-matrix  method.  This 
is  probably  due  to  two  things: 

1.  The  ring-shell  transducers  are  operated  at  resonance  (controlled),  where  the  spherical 
shells  were  not. 

2.  There  is  no  evidence  of  a  spurious  (uncontrolled)  resonance  at  the  operating  frequency 
for  the  ring-shell  transducers,  but  there  was  for  the  spherical  shell  transducers. 


Seven-Element  Array  Configuration 

Figure  5.15  depicts  the  3D  normalized  far-field  beam  pattern  of  a  linear  array  with  seven  evenly- 
spaced,  identical,  ring-shell  transducer  elements  at  “ka=\”  (/= 605  Hz).  The  inter-element  spac¬ 
ing  is  again  equal  to  A/2.  The  width  of  the  main  lobe  is  reduced,  and  the  directivity  is  enhanced 
in  this  configuration  compared  to  the  three  element  array.  Three  conical  nodal  surfaces  can 
clearly  be  seen  in  the  figure. 

Figure  5.16  is  the  normalized  beam  pattern  of  the  seven  evenly-spaced,  identical,  ring-shell 
transducers.  The  inter-element  spacing  is  equal  to  A/2.  Figures  5.16b  and  5.16c  compare  the 
patterns  obtained  from  the  T-matrix  method  versus  the  point  source  calculation.  Although  the 
beam  patterns  of  the  T-matrix  seem  almost  the  same  as  the  beam  patterns  obtained  using  point 
sources  near  broadside,  Figure  5.16c  shows  beam  patterns  obtained  using  T-matrices  do  not 
produce  nulls  as  do  the  point  sources.  Also,  a  2  dB  reduction  was  observed  at  end  fire  for  the 
T-matrix  compared  to  point  sources. 
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Eight-Element  Array  Configuration 

Figure  5.17  depicts  the  3D  normalized  far-held  beam  pattern  of  the  linear  array  with  eight 
evenly-spaced,  identical,  ring-shell  transducer  elements  at  “ka=  1”  (/= 605  Hz).  The  inter¬ 
element  spacing  is  equal  to  A/2.  Four  conical  nodal  surfaces  are  clearly  seen  from  the  figure. 

Although  the  difference  between  the  two  methods  is  practically  indistinguishable  at  broadside, 
Figure  5.18  shows  an  interaction  effect  near  the  end  fire  in  the  detailed  plot  of  amplitudes, 
shown  in  Figure  5.18c,  demonstrating  the  T-matrix  method  does  not  produce  nulls  at  end  fire 
while  the  point  source  method  does.  Also  it  should  be  noted  the  adjacent  nulls  do  not  exist  for 
the  T-matrix  method. 

Three-Element  Array  Configuration  With  Different  Inter-Element  Spacing  Distance 

Figure  5.19  illustrates  the  effect  of  an  inter-element  spacing  on  the  three-element  array  with  “ka” 
equal  to  one.  The  inter-element  spacing  distances  are  arbitrarily  selected  between  A/2-A/4. 
Figures  5.19a  through  5.19i  show  that  decreasing  the  inter-element  distances  increases  array 
interaction.  The  largest  difference  is  observed  for  the  configuration  of  the  A/3  (Figure  5.19f). 
The  difference  between  the  product  theory  and  the  T-matrix  method  are  clearly  observed  near 
the  end  fire  of  this  configuration. 

Array  Configurations  With  Different  Numbers  of  Elements 

To  determine  the  effect  of  the  number  of  elements  in  an  array,  different  configurations  were 
examined.  Figures  5.20  and  5.21  represent  the  far-field  beam  patterns  for  different  numbers  N 
of  elements.  The  elements  are  spaced  equally,  and  inter-element  spacing  is  fixed  at  A/2.  For 
convenience,  even  numbers  of  elements  and  odd  numbers  of  elements  are  displayed  in  separate 
figures. 

Figure  5.21  shows  the  configurations  for  the  even  numbers.  At  end  fire,  an  adjacent  null  is 
presented  with  an  even  number  of  elements  using  the  point  sources  due  to  the  axis  of  symmetry; 
while  the  beam  pattern  obtained  using  the  T-matrix  shows  array  interaction  effects  and  produces 
a  beam  pattern  that  does  not  terminate  at  adjacent  nulls  at  end  fire. 

Figures  5.22  and  5.23  demonstrate  the  3D  normalized  far-field  beam  pattern  for  the  different 
configurations  corresponding  to  figure  Figure  5.20  and  Figure  5.21.  The  directionality  of  the 
array  increases  with  N.  Almost  all  of  the  estimated  conical  nodal  surfaces  by  point  sources  are 
clearly  observed. 
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Summary  of  Results  for  Linear  Arrays  of  Ring-Shell  Transducers 

1.  The  figures  illustrate  the  difference  between  the  two  methods  is  not  remarkable  except 
for  the  end  fire  direction.  Unlike  the  linear  array  of  spherical  thin-shell  transducer,  the 
‘shading”  effect  disappears  for  this  configuration.  Only  the  shadowing  effect  is  produced 
near  the  end  fire,  and  it  is  remarkable  for  even  numbered  configurations. 

2.  The  most  significant  scattering  effect  was  observed  when  the  inter-element  spacing  was 
A/3.  A  spacing  smaller  than  A/3  resulted  in  local  maxima  on  axis  where  the  spherical 
shell  and  point  source  arrays  both  show  local  minima.  The  relative  consistency  of  axial 
to  broadside  amplitude  for  the  ring-shell  array  is  probably  due  to  its  being  generated  at 
resonance. 

3.  The  Figures  5.8  and  5.19  showed  that  for  the  spherical  thin-shell  array,  the  qualitative 
shape  of  the  beam  patterns  is  the  same,  e.g.,  the  angles  of  local  maxima  and  minima. 
However,  for  the  ring-shell  array,  regardless  of  the  spacing,  the  on-axis  response  appears 
to  be  consistently  about  7  dB  smaller  than  the  broadside  response. 

4.  The  least  desirable  array  configuration  of  the  ring-shells  was  observed  for  A/3  inter 
element- spacing  distance  (Figure  5.19f). 
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3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=3  and  d  =>/2 


Figure  5.13:  3D  front-left  view  of  normalized  far  field  beam  pattern  of  vertical  linear  array  with 
three  ring-shell  transducer  elements  obtained  with  the  T-matrix  method, inter-element  spacing 
d=X/ 2,  /= 605  Hz  and  ka=  1 . 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=W2 


(a)  (b)  (c) 


Figure  5.14:  Normalized  far-field  vertical  beam  pattern  for  the  modeled  linear  array  with  three 
identical  ring-shell  transducer  elements  lying  along  polar  axis  (6  =  0),  inter-element  spacing 
d=X /2,  /= 605  Hz  and  ka=  1 .  (a)  Orientation  of  the  ring-shell  transducers  (b)  A  polar  plot  of 
the  normalized  magnitudes  of  the  vertical  far-field  beam  pattern  (c)  A  plot  of  the  normalized 
magnitudes  of  the  vertical  far-field  beam  pattern 
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3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=7  and  d  -V2 


Figure  5.15:  3D  front-left  view  of  normalized  far  field  beam  pattern  of  vertical  linear  array  with 
seven  ring-shell  transducer  elements  obtained  with  the  T-matrix  method,  inter-element  spacing 
d=X/ 2,  /= 605  Hz  and  ka=  1. 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=7  d=M2 


(a)  (b)  (c) 


Figure  5.16:  Normalized  far-field  vertical  beam  pattern  for  the  modeled  linear  array  with  seven 
identical  ring-shell  transducer  elements  lying  along  polar  axis  (6  =  0),  inter-element  spacing 
d=X/2,  f= 605  Hz  and  ka=  1.  (a)  Orientation  of  the  ring-shell  transducers  (b)  A  polar  plot  of 
the  normalized  magnitudes  of  the  vertical  far-field  beam  pattern  (c)  A  plot  of  the  normalized 
magnitudes  of  the  vertical  far-field  beam  pattern 
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3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=8  and  d  =>/2 


Figure  5.17:  3D  front-left  view  of  normalized  far  field  beam  pattern  of  vertical  linear  array  with 
eight  ring-shell  transducer  elements  obtained  with  the  T-matrix  method,  inter-element  spacing 
d=X/ 2,  /= 605  Hz  and  ka=  1. 


(a) 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=8  d=M2 


(b) 


(c) 


Figure  5.18:  Normalized  far-field  vertical  beam  pattern  for  the  modeled  linear  array  with  eight 
identical  ring-shell  transducer  elements  lying  along  polar  axis  (6  =  0),  inter-element  spacing 
d=X/2,  f= 605  Hz  and  ka=  1.  (a)  Orientation  of  the  ring-shell  transducers  (b)  A  polar  plot  of 
the  normalized  magnitudes  of  the  vertical  far-field  beam  pattern  (c)  A  plot  of  the  normalized 
magnitudes  of  the  vertical  far-field  beam  pattern 
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NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/2.2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/2.4 


(a)  d=X/2 


(b)  d=X/  2.2 


(c)  d=X/ 2.4 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=)/2.e  NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/2.8  NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=W3 


(d)  d=X/ 2.6 


(e)  d=X/ 2.8 


(f)  d=X/  3 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/3.2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/3.4  NORM  FF  VERT.  BEAM  PATTERN  in  dB  N=3  d=>/4 


(g)  d=X  /3 .2  (h)  d=X /3 .4  (i)  d=X /4 

Figure  5.19:  The  comparison  of  the  far-held  vertical  beam  patterns  of  a  modeled  evenly  spaced 
three-element  linear  array  of  the  ring-shell  at  /= 605  Hz  and  ka= 1  with  different  inter-element 
spacing  distance  (a-i)  Comparison  of  polar  plot  of  the  normalized  vertical  far-held  beam  pattern 
of  T-matrix  method  against  polar  plot  of  the  normalized  vertical  far-held  beam  pattern  of  the 
point  sources. 
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NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=3  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=5  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=7  and  d=>/2 


(a)  N= 3 


(b)  N= 5 


(c)  N=1 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=9  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=11  and  d=M2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=13  and  d=M2 


(d)  N= 9 


(e)  N=  11 


(f)  N=13 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=15  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=17  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=19  and  d=>/2 


(g)N=15 


(h)  N=17 


(i)  N=19 


Figure  5.20:  The  vertical  far-field  beam  patterns  for  a  modeled  vertical  linear  array  composed  of 
an  odd  number  N  of  evenly  spaced  flextensional  transducer  elements  lying  along  the  polar  axis 
(6  =  0)  at  A/2,  f= 605  Hz  and  ka=  1  (a-i)  Comparison  of  cartesian  plot  of  the  normalized  vertical 
far-field  beam  pattern  of  T-matrix  method  against  cartesian  plot  of  the  normalized  vertical  far- 
field  beam  pattern  of  the  point  sources. 
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NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=4  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=6  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=8  and  d=>/2 


(a)  N= 4 


(b)  N= 6 


(c)  iV=8 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=10  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=12  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=14  and  d=>/2 


(d)  N=  10 


(e)  N=  12 


(f)  N=  14 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=16  and  d=>/2  NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=18  and  d=>/2 


NORM  FF  VERT.  BEAM  PATTERN  in  dB  on  Z  AXiS  N=20  and  d=>/2 


(g)  N=16 


(h)  iV=18 


(i)  N=20 


Figure  5.21:  The  vertical  far-field  beam  patterns  for  a  modeled  vertical  linear  array  composed 
of  an  even  number  N  of  evenly  spaced  spherical  thin-shell  transducer  elements  lying  along  the 
polar  axis  (0  =  0)  at  A/2,  /= 605  Hz  and  ka=  1  (a-i)  Comparison  of  cartesian  plot  of  the  normal¬ 
ized  vertical  far-field  beam  pattern  of  T-matrix  method  against  cartesian  plot  of  the  normalized 
vertical  far-field  beam  pattern  of  the  point  sources. 
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3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=3  and  d  =>/2 


(a)  N= 3 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=9  and  d  =>/2 


(d)  N= 9 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=15  and  d  =V2 


(g)  A/=15 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=5  and  d  =>/2 


(b)  N= 5 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=11  and  d  =V2 


(e)  N=  11 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=17  and  d  =V2 


(h)  N=\l 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=7  and  d  =>/2 


(c)  N=1 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=13  and  d  =V2 


(f)  N=  13 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=19  and  d  =W2 


(i)  Af=19 


Figure  5.22:  The  3D  normalized  far-held  beam  patterns  for  a  modeled  vertical  linear  array 
composed  of  an  odd  number  N  of  evenly  spaced  spherical  thin-shell  transducer  elements  lying 
along  the  polar  axis  (0  =  0)  at  A/2,  /= 605  Hz  and  ka=  1  (a-i)  3D  front-left  view  of  normalized 
far-held  beam  pattern  of  vertical  linear  array 
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(h)  7V=  1 8 

(g)N=16  qualitative 


3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=1,  N=20  and  d  =>/2 


(i)  iV=20 


Figure  5.23:  The  3D  normalized  far-held  beam  patterns  for  a  modeled  vertical  linear  array 
composed  of  an  even  number  N  of  evenly  spaced  spherical  thin-shell  transducer  elements  lying 
along  the  polar  axis  (6  =  0)  at  A/2,  /= 605  Hz  and  ka=  1  (a-i)  3D  front-left  view  of  normalized 
far-held  beam  pattern  of  vertical  linear  array. 
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5.7  General  Conclusion 


The  characteristics  of  beam  patterns  were  obtained  using  the  T-matrix  method  for  modeling  a 
vertical  linear  array  comprised  of  spherical  thin  shell  transducer  elements  and  a  vertical  linear 
array  comprised  of  ring-shell  transducer  elements  at  “ka= 1.”  Results  obtained  with  the  T-matrix 
were  compared  with  for  point  source  arrays.  While,  both  approaches  are  in  agreement  near 
broadside,  they  differ  from  each  other  near  and  at  the  end  fire  directions.  This  is  expected 
because  the  point  sources  ignore  array  element  interaction  while  the  T-matrix  method  considers 
this  affect.  Such  interactions,  if  they  stimulate  uncontrolled  resonances,  as  is  the  case  for  the 
spherical  shell  transducer,  can  result  in  large  deviations  from  point-source  behavior. 

Array  element  interactions  cause  the  degradation  or  the  “shadowing”  in  amplitude  of  beam 
pattern  amplitudes  near  end  fire.  The  shadowing  or  degradation  is  induced  by  scattering  since 
the  scattering  produces  a  destructive  interference  in  an  array  and  changes  the  radiation  away 
from  broadside.  As  a  result  of  these  interactions,  side  lobes  smear  the  minor  lobes  and  shadow 
the  adjacent  null.  In  particular,  this  effect  can  be  seen  in  3D  beam  patterns.  The  nodal  surfaces 
disappear  or  interfere  with  the  closest  side  lobes,  resulting  in  the  expected  nodal  surface  not 
terminating  at  an  adjacent  null  (e.g..  Figure  5.4,  Figure  5.13). 

To  investigate  the  scattering  characteristics,  the  pressure  amplitudes  (e.g,  outgoing,  scattered, 
radiated)  are  normalized  with  the  monopole  radiated  pressure  amplitude,  P(}0.  The  normalized 
pressure  amplitudes  are  multiplied  by  phase  factors  of  the  far-field  Hankel  function,  i^n+l\  in 
order  to  see  the  interaction  effect  on  the  radiation  pattern  explicitly  (equation  (5.5)).  The  nor¬ 
malized  pressure  coefficients  are  listed  in  Tables  5.3  through  5.8.  The  scattering  characteristics 
can  be  deduced  from  these  tables.  The  linear  quadrupole  component  dominates  the  scattering  in 
the  array  consisting  of  the  spherical  shells  while  monopole  scattering  dominates  the  other  com¬ 
ponents  of  the  ring-shell  array.  This  is  entirely  consistent  with  the  previously  mentioned  facts 
that  the  ring-shell  transducers  are  operated  at  (mostly  monopole)  resonance,  where  the  spher¬ 
ical  shell  transducers  are  not,  and  there  exists  an  uncontrolled  linear  quadrupole  resonance  of 
the  spherical  shell  transducers  at  almost  exactly  the  chosen  operating  frequency.  The  differing 
scattering  characteristics  demonstrate  themselves  in  Figures  5.5  and  5.16.  The  increasing  side 
lobes  are  observed  at  end  fire  of  the  spherical  thin-shells,  while  the  side  lobes  of  the  ring-shell 
degrade. 

For  further  investigation  to  understand  the  factors  affecting  the  array  interaction,  the  inter¬ 
element  spacing  distance  was  allowed  to  vary  between  A/2  — A/4.  Resultant  figures  show 


89 


the  effect  of  inter-element  spacing  distances  on  array  sound  pressure  levels.  In  particular,  when 
the  inter-element  spacing  distance  is  equal  to  the  ratio  of  wavelength  to  the  number  of  elements, 
X/N,  significant  changes  are  observed  in  the  beam  patterns  of  Figures  5.19  and  5.8.  The  ar¬ 
ray  interaction  can  be  controlled  by  adjusting  the  inter-element  spacing  distance  in  order  to 
minimize  this  interaction 


The  second  factor  considered  was  the  affect  the  number  of  elements  has  on  beam  pattern  for 
an  array.  This  array  interaction  effect  is  attenuated  as  the  number  of  elements  is  increased.  In 
conclusion,  an  optimal  number  of  elements  can  be  determined  by  analyses  of  beam  patterns  to 
enhance  or  reduce  array  element  interactions. 


N=3 

ka=l  f=474  Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 

0,0 

0.978 

1.000 

0.027 

0.86 

0.00 

146.65 

1 

1,0 

0.068 

0.000 

0.068 

122.39 

90.00 

122.39 

2,0 

1.178 

0.000 

1.178 

320.33 

180.00 

320.33 

0,0 

0.968 

1.000 

0.069 

3.55 

0.00 

119.79 

9 

1,0 

0.000 

0.000 

0.000 

-64.99 

90.00 

-64.99 

2,0 

1.214 

0.000 

1.214 

358.34 

180.00 

358.34 

0,0 

0.978 

1.000 

0.027 

0.86 

0.00 

146.65 

9 

1,0 

0.068 

0.000 

0.068 

-57.61 

90.00 

-57.61 

2,0 

1.178 

0.000 

1.178 

320.33 

180.00 

320.33 

Table  5.3:  The  normalized  pressure  coefficient  of  the  linear  array  consisting  of  spherical  three 
thin-shell  transducers 


N=3 

ka=l  f  =60 5  Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 

0,0 

0.855 

1.000 

0.159 

4.05 

0.00 

157.66 

1 

1,0 

0.031 

0.000 

0.031 

10.46 

90.00 

10.46 

2,0 

0.090 

0.104 

0.014 

176.26 

174.59 

343.96 

0,0 

0.918 

1.000 

0.466 

27.73 

0.00 

113.70 

9 

1,0 

0.000 

0.000 

0.000 

-76.86 

90.00 

-76.86 

z 

2,0 

0.090 

0.104 

0.039 

196.55 

174.59 

296.36 

0,0 

0.855 

1.000 

0.159 

4.06 

0.00 

157.60 

9 

1,0 

0.031 

0.000 

0.031 

190.70 

90.00 

190.70 

2,0 

0.090 

0.104 

0.014 

176.27 

174.59 

343.86 

Table  5.4:  The  normalized  pressure  coefficient  of  the  linear  array  consisting  of  three  ring-shell 
transducers 
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N=4 

ka=lf=474Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 

0,0 

0.984 

1.000 

0.031 

1.54 

0.00 

121.45 

1 

1,0 

0.080 

0.000 

0.080 

104.79 

90.00 

104.79 

2,0 

0.885 

0.000 

0.885 

315.54 

180.00 

315.54 

0,0 

0.970 

1.000 

0.054 

2.65 

0.00 

124.43 

9 

1,0 

0.029 

0.000 

0.029 

18.92 

90.00 

18.92 

2,0 

1.249 

0.000 

1.249 

336.80 

180.00 

336.80 

0,0 

0.970 

1.000 

0.054 

2.65 

0.00 

124.43 

1,0 

0.029 

0.000 

0.029 

198.92 

90.00 

198.92 

2,0 

1.249 

0.000 

1.249 

336.80 

180.00 

336.80 

0,0 

0.984 

1.000 

0.031 

1.54 

0.00 

121.45 

4 

1,0 

0.080 

0.000 

0.080 

-75.21 

90.00 

-75.21 

2,0 

0.885 

0.000 

0.885 

315.54 

180.00 

315.54 

Table  5.5:  The  normalized  pressure  coefficient  of  the  linear  array  consisting  of  four  spherical 
thin-shell  transducers 


N=4 

ka=lf=605  Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 

0,0 

0.913 

1.000 

0.196 

10.53 

0.00 

121.42 

1 

1,0 

0.042 

0.000 

0.042 

-17.41 

90.00 

-17.41 

2,0 

0.094 

0.104 

0.016 

182.29 

174.59 

303.89 

0,0 

0.828 

1.000 

0.349 

19.23 

0.00 

128.73 

9 

1,0 

0.032 

0.000 

0.032 

146.43 

90.00 

146.43 

z 

2,0 

0.084 

0.104 

0.030 

187.95 

174.59 

313.81 

0,0 

0.827 

1.000 

0.349 

19.23 

0.00 

128.77 

1,0 

0.032 

0.000 

0.032 

-33.90 

90.00 

-33.90 

2,0 

0.084 

0.104 

0.030 

187.95 

174.59 

313.87 

0,0 

0.914 

1.000 

0.196 

10.54 

0.00 

121.35 

4 

1,0 

0.042 

0.000 

0.042 

162.81 

90.00 

162.81 

■ 

2,0 

0.094 

0.104 

0.016 

182.30 

174.59 

303.76 

Table  5.6:  The  normalized  pressure  coefficient  of  the  linear  array  consisting  of  four  ring-shell 
transducers 
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N=8 

ka=l  f  =474  Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 

0,0 

0.980 

1.000 

0.032 

1.44 

0.00 

129.79 

1 

1,0 

0.084 

0.000 

0.084 

112.28 

90.00 

112.28 

2,0 

1.033 

0.000 

1.033 

324.51 

180.00 

324.51 

0,0 

0.968 

1.000 

0.058 

2.83 

0.00 

125.04 

9 

1,0 

0.027 

0.000 

0.027 

-22.51 

90.00 

-22.51 

2,0 

1.259 

0.000 

1.259 

345.11 

180.00 

345.11 

0,0 

0.974 

1.000 

0.045 

2.15 

0.00 

126.20 

1,0 

0.033 

0.000 

0.033 

254.49 

90.00 

254.49 

2,0 

1.118 

0.000 

1.118 

329.58 

180.00 

329.58 

0,0 

0.985 

1.000 

0.039 

2.08 

0.00 

113.39 

4 

1,0 

0.051 

0.000 

0.051 

-75.03 

90.00 

-75.03 

2,0 

0.769 

0.000 

0.769 

305.41 

180.00 

305.41 

0,0 

0.985 

1.000 

0.039 

2.08 

0.00 

113.39 

1,0 

0.051 

0.000 

0.051 

104.97 

90.00 

104.97 

w/ 

2,0 

0.769 

0.000 

0.769 

305.41 

180.00 

305.41 

0,0 

0.974 

1.000 

0.045 

2.15 

0.00 

126.20 

1,0 

0.033 

0.000 

0.033 

74.49 

90.00 

74.49 

u 

2,0 

1.118 

0.000 

1.118 

329.58 

180.00 

329.58 

0,0 

0.968 

1.000 

0.058 

2.83 

0.00 

125.04 

7 

1,0 

0.027 

0.000 

0.027 

157.49 

90.00 

157.49 

2,0 

1.259 

0.000 

1.259 

345.11 

180.00 

345.11 

0,0 

0.980 

1.000 

0.032 

1.44 

0.00 

129.79 

8 

1,0 

0.084 

0.000 

0.084 

-67.72 

90.00 

-67.72 

2,0 

1.033 

0.000 

1.033 

324.51 

180.00 

324.51 

Table  5.7:  The  normalized  pressure  coefficient  of  the  linear  array  consisting  of  eight  spherical 
thin-shell  transducers 
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N=8 

ka=l  f  =605  Hz 

Magnitude 

Phase 

(degree) 

Nu. 

Comp. 

Outgoing 

Radiated 

Scattered 

Outgoing 

Radiated 

Scattered 

0,0 

0.902 

1.000 

0.178 

8.99 

0.00 

127.67 

1 

1,0 

0.039 

0.000 

0.039 

-13.03 

90.00 

-13.03 

2,0 

0.094 

0.104 

0.015 

180.85 

174.59 

311.18 

0,0 

0.846 

1.000 

0.380 

21.78 

0.00 

124.38 

9 

1,0 

0.024 

0.000 

0.024 

145.83 

90.00 

145.83 

2,0 

0.085 

0.104 

0.032 

190.45 

174.59 

308.74 

0,0 

0.824 

1.000 

0.290 

14.62 

0.00 

134.21 

1,0 

0.020 

0.000 

0.020 

-48.27 

90.00 

-48.27 

2,0 

0.085 

0.104 

0.024 

183.97 

174.59 

320.30 

0,0 

0.841 

1.000 

0.324 

17.69 

0.00 

127.81 

4 

1,0 

0.017 

0.000 

0.017 

125.19 

90.00 

125.19 

2,0 

0.085 

0.104 

0.027 

186.88 

174.59 

312.79 

0,0 

0.841 

1.000 

0.324 

17.69 

0.00 

127.82 

1,0 

0.018 

0.000 

0.018 

-54.94 

90.00 

-54.94 

w/ 

2,0 

0.085 

0.104 

0.027 

186.88 

174.59 

312.81 

0,0 

0.824 

1.000 

0.290 

14.62 

0.00 

134.19 

1,0 

0.019 

0.000 

0.019 

131.84 

90.00 

131.84 

u 

2,0 

0.085 

0.104 

0.024 

183.97 

174.59 

320.27 

0,0 

0.846 

1.000 

0.380 

21.78 

0.00 

124.40 

7 

1,0 

0.025 

0.000 

0.025 

-34.68 

90.00 

-34.68 

2,0 

0.085 

0.104 

0.032 

190.45 

174.59 

308.77 

0,0 

0.902 

1.000 

0.178 

9.00 

0.00 

127.60 

8 

1,0 

0.039 

0.000 

0.039 

167.19 

90.00 

167.19 

2,0 

0.094 

0.104 

0.015 

180.86 

174.59 

311.05 

Table  5.8:  The  normalized  pressure  coefficient  of  the  linear  array  consisting  of  eight  ring-shell 
transducers 
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CHAPTER  6: 

Conclusion  and  Future  Work 


6.1  Conclusion 

This  thesis  is  the  first  successful  practical  demonstration  of  the  T-matrix  method  applied  to  an 
active  sonar  array  for  which  a  finite-element  model  was  employed  to  compute  the  scattering 
properties  of  a  single  transducer. 

A  piezoelectric  thin  spherical  shell  transducer  and  a  class  V  flextensional  transducer  were  mod¬ 
eled  using  the  COMSOL  Multiphysics  computer  code.  Linear  array  models  of  the  developed 
transducers  were  created  using  the  T-matrix  method. 

Firstly,  we  modeled  a  radially  polarized  piezoelectric  spherical  thin-shell  transducer.  Validation 
of  the  model  was  done  by  comparing  the  computed  values  of  the  free-field  radiated  pressure 
amplitude  against  those  of  an  approximate  analytical  model  derived  from  an  equivalent  circuit 
model  of  the  spherical  transducer.  The  results  agreed  within  3.7  percent,  which  is  acceptable 
for  validation  of  the  finite  element  code. 

As  a  second  transducer  type,  a  class  V  flextensional  ring-shell  transducer  was  selected  in  part 
because  data  are  available  from  a  bi-planar  billboard  array  tested  in  the  1990s  [1].  The  model 
of  the  ring-shell  flextensional  transducer  was  developed  in  COMSOL  Multiphysics.  Computed 
radiation  characteristics  are  in  excellent  agreement  with  the  manufacturer’s  data.  The  primary 
resonance  frequency  computed  is  only  0.8  percent  lower  than  that  of  the  manufacturer’s  data, 
and  the  second  resonance  computed  is  only  0.7  percent  lower.  The  calculated  source  level  (SL) 
for  a  3000  Vrms  drive  for  the  primary  resonance  frequency  was  found  to  be  1.4  dB  lower  than 
the  manufacturer’s  data  at  the  primary  resonance  frequency  and  3  dB  higher  than  the  manufac¬ 
turer’s  data  at  the  secondary  resonance  frequency.  The  bandwidth  and  the  mechanical  quality 
factor  also  show  excellent  agreement  with  the  manufacturer’s  data. 

After  validation  of  the  model,  the  scattering  T-matrix  and  radiation  coefficients  of  the  developed 
model  were  computed  through  quadrupolar  order  using  COMSOL  for  “ ka ”  values  equal  to  one. 
This  same  value  was  used  in  the  T-matrix  method  to  model  the  linear  array. 

Two  linear  array  models  were  created.  Beam  patterns  were  computed  for  different  configura- 
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tions.  The  computed  beam  pattern  results  were  compared  against  those  using  point  source  array 
theory. 

The  most  significant  differences  were  observed  between  the  two  solutions  at  end  fire.  This  is 
expected  because  the  T-matrix  method  incorporates  the  array  interaction  while  the  point  source 
array  theory  does  not.  In  particular,  we  observed  that  the  nulls  at  end  fire  are  not  completed, 
and  they  are  shadowed  by  scattering.  The  T-  matrix  method  accounts  for  the  inter-element  array 
interaction  and  simulates  the  array  performance  more  realistically  than  the  point  source  array 
theory. 

The  T-matrix  method  can  be  used  as  a  versatile  tool  in  the  design  phase  of  an  array  with  high- 
power  active  transducer  elements.  That  is  because  the  multi-interaction  is  accounted  for  in  the 
methodology.  The  accuracy  of  the  T-matrix  method  depends  upon  the  finite-element  model  of 
the  single  transducer  used  to  create  the  array.  COMSOL  Multiphysics  produces  a  good  model 
for  a  single  transducer,  thereby  computing  the  scattering  and  radiation  characteristics  of  arrays 
comprised  of  these  transducers. 

It  should  be  noted,  the  developed  finite  element  model  was  built  for  the  common  case  ka= 1. 
For  higher  ka  values  refinements  should  be  made  in  the  finite  element  model.  Possible  im¬ 
provements  were  investigated  for  the  ring-shell  transducer  in  Chapter  4.4.  When  the  Perfectly 
Matched  Layer  (PML)  boundary  condition  was  applied,  the  model  performed  better  at  frequen¬ 
cies  higher  than  1600  Hz  in  determining  the  SL. 

6.1.1  Future  Work 

As  a  result  of  the  research  done  on  the  T-matrix  method  application  for  array  modeling,  several 
other  interesting  topics  for  study  have  arisen. 

In  comparison  to  an  approximate  method,  it  is  relatively  easy  to  employ  the  T-matrix  method 
in  modeling  active  sonar  array  element  interaction,  and  so  this  method  should  become  the  pre¬ 
ferred  method  for  active  sonar  array  performance  prediction.  However,  more  practical  examples 
are  needed  to  demonstrate  its  applicability  and  superiority  for  an  actual  densely  packed,  high 
power  active  sonar  array. 

In  order  to  test  the  accuracy  of  the  T-matrix  method,  we  are  encouraged  to  develop  an  array 
model  of  the  bi-planar  array  composed  of  24  Sparton  of  Canada  600  Hz  class  flextensional 
ring-shell  transducers  that  were  tested  in  the  1990s  and  compare  the  data  to  our  computed 
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results  using  the  T-matrix  method. 

As  an  alternative  to  computation  of  the  mutual  acoustic  radiation  with  Pritchard’s  approxima¬ 
tion,  the  development  of  a  mutual  acoustic  radiation  impedance  using  a  T-matrix  description 
can  be  carried  out  by  future  research  in  order  to  evaluate  a  projector’s  internal  stress. 

The  accuracy  of  PML  applied  in  the  3D  finite  element  model  can  also  be  investigated  for  im¬ 
proving  the  performance  of  the  3D  model  at  high  frequencies. 
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APPENDIX  A: 

Transforming  the  Elasto-Piezoelectric  Property  Matrix  to 
Model  a  Circumferentially-Polarized  Ring  as 
Radially-Polarized 


The  elasto-piezoelectric  property  matrix  for  a  piezoelectric  ceramic,  taking  stress  T  and  electric 
field  E  as  independent,  is 
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In  a  circumferentially-polarized  ring  for  a  right-handed  coordinate  basis  set,  z  =  1,  r  =  2,  and 
0=3  (axial, radial, circumferential,  respectively).  Then,  denoting  /  for  longitudinal,  s'  for  shear, 
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=  Dr 

0 

0 

0 

^15 

0 

0 

0 

pT 

fcll 

0 

E2 

=  Er 

d3 

=  Dp, 

.  ^31 

d3 1 

d33 

0 

0 

0 

0 

0 

PT 

fc33  J 

.  £3 

=  Ep,  _ 

(A.2) 


Now,  suppose  we  want  to  model  this  ring  as  radially  polarized,  while  trying  to  keep  its  elastic 
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and  dielectric  matrix  as  follows. 


Exchange  the  last  two  rows  and  columns  (i.e.,  8  and  9) 


Si 

=  Siz 

r 

4n 

S12 

413 

0 

0 

0 

0 

^31 

0 

'  Ti 

=  Th  ' 

S2 

=  Sir 

sE 

412 

SE 

411 

413 

0 

0 

0 

0 

^31 

0 

T2 

=  Th 

S3 

=  Si4 

4 13 

413 

,E 

433 

0 

0 

0 

0 

d33 

0 

73 

=  7/0 

s4 

=  Ssz 

0 

0 

0 

,E 

4  44 

0 

0 

0 

0 

d\5 

t4 

=  Esz 

S5 

=  Ssr 

= 

0 

0 

0 

0 

,£ 

4  44 

0 

dis 

0 

0 

t5 

=  Tsr 

S6 

=  SS(j) 

0 

0 

0 

0 

0 

*E 

s66 

0 

0 

0 

t6 

=  Ts(p 

Di 

=  DZ 

0 

0 

0 

0 

dl5 

0 

pT 

fcll 

0 

0 

Ei 

=  EZ 

D2 

=  Dr 

^31 

^31 

d33 

0 

0 

0 

0 

pT 

fc33 

0 

E2 

=  Er 

D3 

=  D(t, 

0 

0 

0 

^15 

0 

0 

0 

0 

£T 
fcll  J 

.  £3 

=  E<t>  . 

(A.3) 


Change  the  correspondence  of  the  matrix  indices  with  coordinate  directions  to  identify  the  new 
polarization  direction  as  the  “3”  direction  (0  =  2,  r  =  3).  Exchanges  rows  and  columns  2  and  3, 
5  and  6,  and  8  and  9 


Si 

=  Siz 

r 

4n 

413 

4 12 

0 

0 

0 

0 

0 

^31 

'  7i 

=  7}z  ' 

s2 

=  S[(f> 

413 

v£ 

433 

413 

0 

0 

0 

0 

0 

d33 

72 

=  7/0 

S3 

=  Sir 

v£ 

412 

v£ 

413 

v£ 

411 

0 

0 

0 

0 

0 

d3\ 

73 

=  Tir 

s4 

=  Ssz 

0 

0 

0 

4  44 

0 

0 

0 

<7l5 

0 

74 

=  7^ 

S5 

=  Ss(f) 

= 

0 

0 

0 

0 

.E 

s66 

0 

0 

0 

0 

7s 

=  7^0 

S6 

=  SSr 

0 

0 

0 

0 

0 

,E 

4  44 

<7l5 

0 

0 

76 

=  7^,- 

Di 

=  E>Z 

0 

0 

0 

0 

0 

<7l5 

£T 

fcll 

0 

0 

El 

=  ^z 

d2 

=  Dft, 

0 

0 

0 

^15 

0 

0 

0 

ft 

fcn 

0 

e2 

=  E<p 

d3 

=  Dr  _ 

.  <731 

<733 

^31 

0 

0 

0 

0 

0 

PT 

fc33  J 

.  £3 

=  Er  _ 

(A.4) 


The  final  form  is  not  expected  to  depend  upon  whether  the  transformed  polarization  direction 
is  chosen,  r  or  z.  To  prove  this,  we  follow  the  same  procedure  as  above,  now  exchanging  </>  for 
the  z  polarization  direction.  We  exchange  rows  and  columns  (7  and  9). 
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Si 

=  Slz 

r 

■'ll 

sE 

412 

413 

0 

0 

0 

^31 

0 

0 

'  Ti 

=  Th  " 

S2 

—  Sir 

sE 

412 

■'ll 

4 13 

0 

0 

0 

d3\ 

0 

0 

Ti 

=  Tir 

S3 

=  Siq 

v£ 

413 

v£ 

413 

v£ 

433 

0 

0 

0 

d33 

0 

0 

T3 

=  Ti<p 

s4 

—  Ssz 

0 

0 

0 

,E 

4  44 

0 

0 

0 

^15 

0 

t4 

=  Esz 

S5 

=  Ssr 

— 

0 

0 

0 

0 

.E 

4  44 

0 

0 

0 

^15 

t5 

=  TSr 

s6 

—  SS(j) 

0 

0 

0 

0 

0 

s66 

0 

0 

0 

t6 

=  Tstf, 

Di 

=  DZ 

^31 

^31 

d33 

0 

0 

0 

?T 

fc33 

0 

0 

Ei 

=  EZ 

d2 

=  D, 

0 

0 

0 

^15 

0 

0 

0 

pT 

fcll 

0 

E2 

=  Er 

D3 

=  Dty 

0 

0 

0 

0 

d\5 

0 

0 

0 

pT 

fcll  J 

.  £3 

=  E<t>  . 

(A.5) 


We  change  the  correspondence  of  matrix  indices  with  coordinate  directions  to  identify  the  new 
polarization  direction  as  the  “3”  direction  (we  now  assigning  (j)  =  1,  z  =  3) 


S3 

=  Siz 

r  v£ 

4n 

v£ 

4 12 

413 

0 

0 

0 

^31 

0 

0 

'  T3 

=  Tlz  ' 

s2 

=  Sir 

4 12 

411 

v£ 

413 

0 

0 

0 

d3l 

0 

0 

t2 

=  Tir 

Si 

=  S14 

v£ 

413 

v£ 

413 

v£ 

433 

0 

0 

0 

d23 

0 

0 

Ti 

=  Tiq 

S4 

=  Ssz 

0 

0 

0 

,E 

4  44 

0 

0 

0 

^15 

0 

T6 

=  Esz 

S5 

=  Ssr 

= 

0 

0 

0 

0 

4  44 

0 

0 

0 

dis 

t5 

=  TSr 

S6 

=  SS(f) 

0 

0 

0 

0 

0 

s66 

0 

0 

0 

t4 

=  Ts(p 

D3 

=  DZ 

<^31 

<^31 

d33 

0 

0 

0 

pT 

b33 

0 

0 

e3 

=  EZ 

d2 

=  Dr 

0 

0 

0 

^15 

0 

0 

0 

pT 

fcll 

0 

e2 

=  Er 

Di 

=  E>(f, 

0 

0 

0 

0 

^15 

0 

0 

0 

pT 

fcll  J 

.  El 

=  E(j,  _ 

(A.6) 


We  exchange  rows  and  columns  1  and  3,  4  and  6,  7  and  9  to  put  them  in  regular  order: 
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Si 

=  Sty 

r 

433 

413 

sE 

413 

0 

0 

0 

0 

0 

d33 

'  TX 

=  Ti^ 

S2 

—  Sir 

413 

sE 

■'ll 

v£ 

4 12 

0 

0 

0 

0 

0 

d3i 

Ti 

=  Tir 

S3 

=  slz 

v£ 

413 

v£ 

412 

v£ 

411 

0 

0 

0 

0 

0 

<*3l 

T3 

=  Th 

s6 

—  Sty 

0 

0 

0 

SE 

466 

0 

0 

0 

0 

0 

t4 

=  Es(j, 

S5 

=  Ssr 

= 

0 

0 

0 

0 

4  44 

0 

<*15 

0 

0 

t5 

=  Tsr 

S4 

=  Ssz 

0 

0 

0 

0 

0 

4  44 

0 

<*15 

0 

t6 

=  Tsz 

Di 

=  Dq 

0 

0 

0 

0 

<*15 

0 

oT 

fcll 

0 

0 

Ei 

=  £> 

d2 

—  Dr 

0 

0 

0 

0 

0 

<*15 

0 

pT 

fcll 

0 

E2 

=  Er 

d3 

=  DZ  . 

.  d33 

<*31 

<*3l 

0 

0 

0 

0 

0 

PT 

fc33  J 

.  £3 

=  EZ  . 

(A.7) 


An  alternative  procedure  to  transform  the  electrical  coupling  is  as  follows.  Instead  of  exchang¬ 
ing  rows  and  columns,  permute  them.  Permuting  the  polarization  direction  to  “r”  (i.e.,  the  “2” 
direction): 


Si 

=  sh  ' 

r 

4n 

v£ 

412 

v£ 

413 

0 

0 

0 

0 

<*31 

0 

'  Ti 

=  Th  ' 

s2 

=  Sir 

v£ 

412 

v£ 

411 

v£ 

413 

0 

0 

0 

0 

<*31 

0 

T2 

=  Tlr 

S3 

=  Slip 

413 

413 

433 

0 

0 

0 

0 

<*33 

0 

T3 

=  Tty 

S4 

=  Ssz 

0 

0 

0 

.E 

4  44 

0 

0 

<*15 

0 

0 

T4 

=  TSZ 

S5 

=  Ssr 

= 

0 

0 

0 

0 

,E 

4  44 

0 

0 

0 

<*15 

t5 

=  Tsr 

S6 

=  Ss<j) 

0 

0 

0 

0 

0 

„E 

466 

0 

0 

0 

t6 

=  Estp 

Di 

=  DZ 

0 

0 

0 

<*15 

0 

0 

pT 

fcll 

0 

0 

El 

=  EZ 

D2 

=  Dr 

<*31 

<*31 

<*33 

0 

0 

0 

0 

pT 

fcll 

0 

El 

=  Er 

D3 

=  E>(j> 

0 

0 

0 

0 

<*15 

0 

0 

0 

pT 

fc33  J 

_e3 

=  E<S>  . 

(A.8) 


Now,  reidentify  coordinate  directions  and  matrix  subscripts: 
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S2 

=  Slz 

r 

■'ll 

412 

413 

0 

0 

0 

0 

^31 

0 

H 

to1 

II 

M 

S3 

—  Sir 

sE 

412 

sE 

■'ll 

4 13 

0 

0 

0 

0 

^31 

0 

II 

"S 

Si 

=  Siq 

v£ 

413 

v£ 

413 

v£ 

433 

0 

0 

0 

0 

d23 

0 

II 

S5 

—  Ssz 

0 

0 

0 

4  44 

0 

0 

d\5 

0 

0 

N 

10 

11 

Se 

=  Ssr 

— 

0 

0 

0 

0 

.E 

4  44 

0 

0 

0 

d\5 

II 

''S 

S4 

Ss(j) 

0 

0 

0 

0 

0 

466 

0 

0 

0 

-e- 

E-? 

II 

D2 

=  DZ 

0 

0 

0 

^15 

0 

0 

oT 

fcll 

0 

0 

to 

II 

M 

D3 

=  D, 

^31 

^31 

d23 

0 

0 

0 

0 

ft 

fcil 

0 

to 

II 

"S 

Di 

=  Dty 

0 

0 

0 

0 

d\s 

0 

0 

0 

PT 

fc33  J 

II 

"O- 

1 _ 

Now,  put  the  matrix  in  regular  order: 


Si 

=  Sty 

433 

,E 

413 

,E 

412 

0 

0 

0 

0 

0 

d33 

'  T\  = 

Tty 

S2 

=  sh 

„E 

413 

,E 

411 

,E 

413 

0 

0 

0 

0 

0 

d3l 

t2  = 

Tiz 

S3 

=  Sir 

*E 

412 

.E 

413 

.E 

411 

0 

0 

0 

0 

0 

d3i 

?3  = 

Tlr 

S4 

=  Sstj) 

0 

0 

0 

SE 

4  66 

0 

0 

0 

0 

0 

74  = 

7j0 

S5 

=  Ssz 

= 

0 

0 

0 

0 

*E 

4  44 

0 

0 

^15 

0 

t5  = 

Tsz 

S6 

=  Ssr 

0 

0 

0 

0 

0 

,E 

444 

d\5 

0 

0 

76  = 

Tsr 

Di 

=  D,)) 

0 

0 

0 

0 

0 

^15 

PT 

fcll 

0 

0 

Ei  = 

E<f, 

d2 

=  DZ 

0 

0 

0 

0 

^15 

0 

0 

PT 

fcll 

0 

e2  = 

Ez 

D3 

=  Dr 

.  d33 

^31 

<^31 

0 

0 

0 

0 

0 

PT 

fc33  J 

.  E3  = 

Er  _ 

Following  the  permutation  procedure,  now  make  “z”  the  polarization  direction: 


(A.9) 


(A.  10) 
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(A.ll) 


Si 

=  Slz 

SE 

*11 

,E 

S12 

,E 

*13 

0 

0 

0 

d3l 

0 

0 

'  T\ 

=  Tlz  ' 

S2 

=  Sir 

,E 

s12 

SE 

*11 

,E 

*13 

0 

0 

0 

d3l 

0 

0 

t2 

=  Tir 

S3 

=  Sty 

*13 

*E 

*13 

*E 

*33 

0 

0 

0 

d33 

0 

0 

t3 

=  Tty 

s4 

=  Ssz 

0 

0 

0 

VE 

*44 

0 

0 

0 

0 

d\s 

t4 

=  Tsz 

S5 

=  Ssr 

— 

0 

0 

0 

0 

VE 

*44 

0 

0 

^15 

0 

t5 

=  Tsr 

s6 

II 

& 

-0- 

0 

0 

0 

0 

0 

<E 

*66 

0 

0 

0 

t6 

II 

-0- 

Di 

=  DZ 

^31 

^31 

d33 

0 

0 

0 

pT 

fc33 

0 

0 

Ei 

=  EZ 

d2 

=  Dr 

0 

0 

0 

0 

d\$ 

0 

0 

pT 

fcll 

0 

e2 

=  Er 

d3 

=  Dq 

0 

0 

0 

^15 

0 

0 

0 

0 

pT 

fcll  J 

e3 

=  Eff, 

Now,  reidentify  coordinate  directions  and  matrix  subscripts: 


s3 

=  sh  ' 

*11 

„E 

*12 

„E 

*13 

0 

0 

0 

6?31 

0 

0 

'  T3 

=  Tiz  ' 

Si 

=  Sir 

*E 

*12 

.E 

■^11 

*E 

*13 

0 

0 

0 

^31 

0 

0 

Ti 

=  T,r 

S2 

=  Sty 

„E 

S13 

,E 

*13 

,E 

s33 

0 

0 

0 

d33 

0 

0 

t2 

=  Tty 

S() 

=  Ssz 

0 

0 

0 

<E 

*44 

0 

0 

0 

0 

^15 

t6 

=  Tsz 

s4 

=  Ssr 

= 

0 

0 

0 

0 

<E 

*44 

0 

0 

^15 

0 

t4 

=  Tsr 

S5 

=  SS(j) 

0 

0 

0 

0 

0 

sE 

*66 

0 

0 

0 

t5 

=  Tty 

D3 

=  DZ 

^31 

^31 

d33 

0 

0 

0 

pT 

fc33 

0 

0 

e3 

=  EZ 

Di 

=  Dr 

0 

0 

0 

0 

^15 

0 

0 

pT 

fcll 

0 

El 

=  Er 

d2 

=  E>(j) 

0 

0 

0 

^15 

0 

0 

0 

0 

PT 

fcll  J 

e2 

=  E(/) 

(A.  12) 
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Now  put  in  regular  order: 


Si 

=  Sir 

511 

CE 

513 

,E 

512 

0 

0 

0 

0 

0 

d3\ 

'  TX 

=  Tir 

s2 

=  Si^ 

*E 

a13 

*E 

a33 

*E 

a13 

0 

0 

0 

0 

0 

d33 

t2 

=  Ti(j) 

S3 

=  S[z 

a12 

a13 

*E 

^ll 

0 

0 

0 

0 

0 

d3i 

t3 

=  Tiz 

S4 

=  Ssr 

0 

0 

0 

<E 

s44 

0 

0 

d\5 

0 

0 

t4 

=  Tsr 

S5 

II 

& 

■0- 

= 

0 

0 

0 

0 

<E 

s66 

0 

0 

0 

0 

t5 

II 

-0- 

S6 

—  Ssz 

0 

0 

0 

0 

0 

*E 

s44 

0 

d\s 

0 

t6 

—  Esz 

Di 

=  Dr 

0 

0 

0 

d\5 

0 

0 

pT 

fcll 

0 

0 

Ei 

=  Er 

d2 

=  Dty 

0 

0 

0 

0 

0 

^15 

0 

pT 

fcll 

0 

e2 

=  E,/, 

d3 

=  Dz  . 

.  d33 

^31 

<*31 

0 

0 

0 

0 

0 

PT 

fc33  J 

.  E3 

=  EZ  _ 

(A.13) 
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APPENDIX  B: 
Main  Code 


%case:  arbitrary  number  of  element  on  the  Z(polar)  axis 
clear  all 
clc 

element=[2]; 

distance=[2]; 

for  numl=l:length(element) 

for  num2=l  :length(distance) 

clearvars  -except  element  distance  numl  num2 

clc 

mode=2;  %  total  number  of  mode 

K2=(mode+l )'  ‘2;  %  monopole, dipole,  quadrupole 

%  Parameter  are  adjusted  for  ka=l 

radius=0.4;  %radius  of  ring-shell  [m] 

c=1490;  %speed  of  sound  [m/s] 

f=605; 

lambda=c/f;  %wavelength  [m] 
k=2*pi/lambda;  %wave  number 
ka=k*radius;  %  ka=1.0 

d=lambda/distance(num2);  %distance  between  the  projector 
kd=k*d; 

%write  transducer  location  for  line  array 

N=element(numl);  %number  of  projector 

%check  the  number  of  element  whether  odd  or  even  number 

if  mod(N,2)==0  %if  it  is  even  the  spacing  along  the  z  axis  would 

p=-N/2:N/2-l; 

z=(p*d)+d/2;  %transducer  located  through  axes(x,y,z) 
else 

p  =-(N-l)/2: 1  :(N-l)/2; 

z=p*d;  %transducer  located  through  axes(x,y,z) 
end 
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%other  axis  would  be  zero  for  linear  array 

x=zeros(l,N); 

y=zeros(l,N); 

position=[x’  y’  z’]; 

%calculate  the  relative  angle  for  each  transducer 

for  varl=l:N 

for  var2=l:N 

delta_x=0; 

delta_y=0; 

delta_z=z(var  1  )-z(var2) ; 

r=sqrt(delta_x.  ~2+delta_y.  ~2+delta_z.  ~2); 

theta(varl,var2)=acos(delta_z/r); 

phi(varl,var2)=atan2(delta_x,delta_y); 

end; 

end; 

%Set  the  diagonal  angles  zeros(i.e  thetall,theta22,theta33...) 
theta(  1  :N+ 1  :N*N)=0; 

%%  Computation  of  translation  Matrix  Gij  for  transducer  elements 

%  the  modes  of  radiated  transducer  called  as  r,m  and  the  modes  of  the  transducer 

%  where  the  radiated  pressure  defined  as  incoming  pressure  reference  with 

%  it’s  origin  called  as  s,u. 

alfa=zeros(K2*N); 

for  sph_m=l:N  %main  ring-shell 

for  sph_r=l:N  %relative  ring-shell 

kr  1 2=abs(sph_m-sph_r)*kd; 

count2=0; 

for  r=0:mode  %2nd  ring-shell  modes  r,m 
for  m=-r:r 

count2=l+count2;  %index  of  Gij  matrix  column 

count  1=0;  %index  of  Gij  matrix  row 

for  s=0:mode  %lth  ring-shell  modes  s,u 

for  u=-s:s 

count  l=count  1+1; 

%define  the  lower  limit  of  summation 
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tO=abs(r-s); 
tl=abs(m-u); 
if  tl>tO 

if  mod((tl-t0),2)==0 

tO=tl; 

else 

tO=tl+l; 

end; 

end; 

count3=l; 
clear  Gtemp3 
for  t=t0:2:s+r 

%Van  Buren’s  paper  used  for  defining  the  a 
%coeffcient  and  named  this  function  as  "a". 
Gtempl=a(s,t,r,u,m)*shank2(t,krl2); 

%spherical  harmonic  function  calculate  sph. 

%harmonic 

M=m-u;  %order  of  Spherical  Harmonic 

Gtemp2=sphericalharmonic(t,M,theta(sph_m,sph_r),phi(sph_m,sph_r)); 
%Multiplication  of  defined  spherical  harmonic, a 
%coeff  and  sph.  hankel. 

Gtemp3(l,count3)=Gtempl*Gtemp2; 

%increase  the  count3  t  create  a  Btemplc 
%vector. 

count3=count3+l ; 
end; 

%the  translation  matrix  which  showed  the  outgoing  wavw 

%from  2  sphere  as  a  incominwave  to  1  sphere. 

if  krl2==0 

G(count  1  ,count2)=0; 

else 

G(countl  ,count2)=sum(Gtemp3); 

end; 

end; 
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end; 

countl=l; 

end; 

end; 

%locate  the  product  matrix(T- matrix  and  Gij)  as  a  submatrix  into  Alfo  matrix 

load(’Tmatrix_ring.mat’)  %T-MAtrix  computed  by  the  help  of  COMSOL 

pos  1  a=(  1  +length(G)  *(sph_r- 1 )) ; 

pos  lb=(length(G)+length(G)*(sph_r- 1 )); 

pos2a=  1  +length(G)  *  (sph_m- 1 ) ; 

pos2b=(length(G)+length(G)*(sph_m-l)); 

if  sph_r-sph_m==0 

alfa(pos2a:pos2b,posla:poslb)=zeros(length(Tmatrix)); 

else 

alfa(pos2a:pos2b,posla:poslb)=Tmatrix*G; 

end; 

end; 

end; 

%%  Calculate  Beam  Pattern 

clearvars  -except  alfa  N  position  element  distance  lambda  d  k  mode  numl  num2 

load(’ Radiated  Pressure  Amplitude_ringlV.mat’)  %Radiated  pressure  computed  by  the  help  of 

COMSOL 

Pr=repmat(P_radiated,N,l);  %repeat  the  vector  for  the  specified  mode  for  each  projector 

%%  DEFINE  SCATTER  VECTOR 

%express  in  matrix  form  A*P_s=B 

%A  :  I- Alfa  and  B:  Free  Radiated  Pressure  multiply  alfa 

A=eye(length(alfa))-alfa; 

B=alfa*Pr; 

%Solution  for  P_s 
P_s=A\B; 

%%  Total  Pressure  Computation 

%Point  position, spherical  coordinate  system  according  to  the  top  transducer 
sample=30; 

phi=linspaee(0, pi, sample) ; %  degree 
%convert  radian 
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theta=linspace(0,  pi,  sample) ; 

%define  the  position  vectors  and  unit  vectors 
position=position’ ; 

%Summation 

for  ind2=l  :length(theta) 

rhat=  [cos(phi).  *  sin(theta(ind2))  ;sin(theta(ind2)).  *  sin(phi)  ;cos(theta(ind2)) .  *ones(  1 , length  (phi))] ; 
for  indl=l:length(phi) 

Ptotal=0; 

count=0; 

for  proj=l:N  %number  of  projector 

Ptempb=0; 
for  n=():mode 
Ptempa=0; 
for  m=-n:n 
count=count+l; 

%Spherical  Harmonic  Function  changed, because  theta  is  not  scalar 
%anymore 

Ptempa=Ptempa+(Pr(count)+P_s(count))*sphericalharmonic(n,m,theta(ind2),phi(indl));  %  to¬ 
tal  Pressure 
end 

Pternpb=Ptempb+Ptempa*(li)~(n+l);  %  total  Pressure 
end 

Ptotal=Ptotal+Ptempb.  *exp(  1  i.  *k.  *(position(  1  ,proj ).  *rhat(  1  ,ind  1 ). . . 

+position(2,proj ).  *rhat(2,ind  1  )+position(3  ,proj ) .  *rhat(3  ,ind  1 ))); 
end 

Ptotall(indl,ind2)=Ptotal; 

end; 

%Find  the  beam  pattern 
end; 

%normalization 

H=Ptotall/abs(max(max(Ptotall)));  %  the  max  function  takes  the  magnitude  of  the  complex 
value 

b=20*logl0(abs(H));  %Beam  Pattern  in  dB 

shift=50; 
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b  =b+shift; 
b(find(b<0))  =  0; 

%%  Plots  of  Results 

x  =  (b.*(cos(phi)’*sin(theta))); 
y  =  (b.*(sin(phi)’*sin(theta))); 
z  =  b.  *repmat(cos(theta), sample, 1); 
x  1  =-(b.  *  (cos(phi)  ’  *sin(theta))) ; 
y  1  =-  (b.  *  (sin(phi)  ’  *  sin(theta) )) ; 
zl  =b.*repmat(cos(theta), sample,  1); 

%For  convenience  in  linear  Plot  convert  the  angle  unit  from  radian  to  degree 
phi_d=phi*  1 80/pi; 
figure  l=figure 
axis  equal 

set(gca,  ’  FontS  ize  ’,12) 
surf(x,y,z,b) 

hold  on 

surf(xl,yl,zl,b) 

axis  square;  camlight  headlight; 

set(gca,  ’  box  ’ ,  ’  off’ ) ; 

set(gca,  ’  xtick’ ,  [] ,  ’  y  tick’ ,  [] ,  ’  ztick’ ,  [] ) ; 

xlabel(”) 

ylabel(”) 

zlabel(’ ARRAY  AXIS’, ’FontSize’,  18) 

title([’3D  Far  Field  Beam  Pattern  of  Linear  Array, ka=l,  N=’,num2str(N)... 

,’  and  d  =  \lambda/’,num2str(distance(num2))],’FontSize’,15) 

screen=get(0,  ’  ScreenSize  ’ ) ; 

end; 

end; 


B.l  Code  of  the  Used  Functions 

The  following  code  was  used  in  the  main  code  they  are  written  by  Joseph  L.  Day  [13] 
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shank2 


function  hn=shank2(n,x) 

%  Computes  the  Spherical  Hankel  function  of  the  second  kind 
%  Written  by  Joseph  L.  Day 

hn=sqrt(pi/(2*x))*besselh(n+.5,2,x); 


B.1.1  jprime 

function  jp=jprime(n,x) 

%  Computes  the  first  derivative  of  the  Spherical  Bessel  function. 
%  Written  by  Joseph  L.  Day 

jp=sqrt(pi/(2*x))*(besselj(n-.5,x)-(n+l)*besselj(n+.5,x)/(x)); 


B.1.2  fact 

function  y=fact(n) 

%  Computes  the  factorial  of  n 
%  Written  by  Joseph  L.  Day 

y=prod(l:n) ; 


B.1.3  be 

function  y=bc(n,m) 

%  Computes  the  binomial  coefficient 
%  Written  by  Joseph  L.  Day 
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if  m<0 

y=0; 

el  seif  n-m<0 

y-0; 

else 

y=prod(l  :n)/(prod(l  :m)*prod(l  :n-m)); 
end 


B.1.4  a 

function  A=a(s,t,r,u,m) 

%  computes  the  function  a(s,t,r,u,m)  of  the  Addition  Theorem 
%  as  written  in  King  and  Van  Buren 

num  1 =(2*s+ 1 )  *  (2*t+ 1 )  *fact(s-u)  *fact(t-m+u)*fact(r+m)  *fact((s+t+r)/2) ; 
deml=fact((r+t-s)/2)*fact((r+s-t)/2)*  fact((s+t-r)/2)*fact(s+t+r+l); 
term  1  =num  1  /  dem  1 ; 

wmin=max([r-s-t,s-r-t-2*m+2*u,t-s-r+2*u])/2; 

wmax=min([s+t-r,r+t-s-2*m+2*u,r+s-t+2*u])/2; 

sum=0; 


for  w=wmin:wmax 

sum=sum+(- 1 ) "  w*bc(s+t-r,(s+t-r)/2+w)*.. . 

bc(t+r-s,(t+r-s)/2+m-u+w)*bc(s+r-t,(s+r-t)/2-u+w); 

end 

A=(i~(s+t-r))*terml*sum; 
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Initial  Distribution  List 


1 .  Defense  Technical  Information  Center 
Ft.  Bel  voir,  Virginia 

2.  Dudley  Knox  Library 
Naval  Postgraduate  School 
Monterey,  California 
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